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ABSTRACT 


First-order solutions indicate tiiat a forced Keplerian trajectory (FKT) 
obtained by thrust-drag cancellation is as fiiel-efficient as a Hohmann transfer. 
Further analysis has shown that the FKT is not Mayer-optimal. Therefore there 
must exist another trajectory that matches or exceeds the efficiency of the 
Hohmann transfer. The application of this result to the fuel-optimal orbit 
maintenance problem implies that periodic reboosts must be more efficient than an 
FKT profile. This research begins with the formulation of an optimal periodic 
cmith)! (OPC) problem to determine the minimum fuel-reboost strategy. The 
problem is muneiically solved by a spectral collocation method. The optimization 
code is further modified to increase accuracy and reduce sensitivity to initial 
guesses. The results of this effort identified a trajectory for a sample satellite that 
was 3.5% more efficient than an ideal impulsive Hohmann transfer over the same 
period of time. From the optimal code, a maximum thruster size is also 
identifiable for a set of initial conditions. The optimal trajectory can save as much 
as 10% of the propellant budget when compared to finite-bum Hohmann transfers. 
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INTRODUCTION 


The boundary of the Earth’s atmosphere extends well into the operating area of low- 
Earth-orbiting satellites. As a vehicle circles the Earth at over seven kilometers per 
second, the few molecules that exist in the upper layers of the atmosphere are continualfy 
striking the vehicle’s surfece. Over time, these collisions decrease the orbital energy of the 
spacecraft, lowering its orbital altitude to an area filled with higher concentrations of 
molecules. This perturbing force is called atmo^heric drag and is a large factor in 
determining the lifetime of any low-Earth-orbiting satellite. While the density of the 
atmosphere at altitudes where these satellites operate, 200 km to 500 km above the 
Earth’s surfece, is low (~10'" kg/m^), the velocity of the spacecraft is very high (~8 km/s). 
These two fectors combine to produce a feirly significant drag force that operates in 
opposition to the velocity vector. A satellite that Ms to counter these effects will 
continue to lose its orbital altitude, spiraling towards the Earth until it bums-up during 
reentry. To avoid this situation, mission lifetimes are extended by maneuvers that reboost 
the satellite to its original altitude periodically to prevent destructive atmospheric heating. 
This reboost maneuver is the largest source of propellant consumption for many low- 
Earth-orbiting vehicles and mission lifetime depends directly on the amount of fuel a 
spacecraft can carry for reboost. This thesis examines this process and attempts to 
maximize the satellite fuel efficiency by determining the optimal reboost trajectory through 
the use of optimal periodic control. 

As a spacecraft circles the Earth, drag continually removes energy from the orbit. This 
causes an instantaneous reduction of the vehicle’s orbital velocity. The change in velocity 
alters the original orbital shape, increasing the eccentricity. In a sense, the satellite enters 
into a “steeper” orbit that in turn increases the orbital velocity. This is known as the drag 
paradox, where the removal of orbital energy causes a satellite to travel fester. Since drag 
is higher at greater densities and velocities, the perturbing force of drag increases. The 
result is that the satellite follows a spiral path, which increases in velocity as the altitude 
decreases. Eventually drag will grow vmtil significant atmospheric heating begins to 
occur. The reboost maneuver avoids a premature end to the satellite’s lifetime by 
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reboosting it to a higher altitude, a region of reduced drag. Figure 1.1 shows this 
trajectory pattern and orbit transfer process. 


DRAG INDUCED SPIRAL AND 
REBOOST MANEUVER 



Figure 1-1 

The most common type of orbital transfer, mainly in terras of mission planning, is the 
Hohmann transfer. This maneuver begins with the satellite in a circular orbit at its lowest 
desired height. To reboost the vehicle, an impulsive bum is made which places the 
satellite into an elliptical transfer orbit. This transfer orbit has two important geometric 
properties. First, the perigee of this orbit is located at the position of the initial impulsive 
bum and secondly, the apogee of the orbit has an altitude that is equal to the desired final 
altitude. When the satellite reaches apogee, another impulsive bum is made to circularize 
the final orbit. 

Previous work has examined the Hohmann transfer against other trajectories [Ref 
1,2]. Using first order estimates of different trajectories, optimal control theory states that 
the Hohmann transfer is not the most fuel-efficient solution [Ref 3,4]. Thus there must 
exist a trajectory that maintains the satellite in a desired orbital band but with less 
consumed propellant. The previous thesis efforts at the Naval Postgraduate School have 



examined several orbital transfer methods but did not surpass the efficiency of the 
Hohmann transfer. This effort takes a more rigorous approach in searching for the 
optimal solution to the orbit maintenance problem Optimal periodic control (OPC) theory 
is applied to the system, which is governed by a set of equations of motion. Numerically 
solving the OPC problem gives a solution and a trajectory that can be adequately 
compared to the Hohmann transfer. 

Instead of guessing a trajectory first and then comparing its performance against the 
Hohmann, this method mathematically derives an optimal trajectory that then can be 
compared. The goals of this thesis are twofold. First, this thesis will analyze the orbital 
transfer by numerically solving the OPC problem Second, the thesis will identify a 
trajectory that is more fuel-efficient than currently used orbital transfer techniques. 


Launch Vehicle 

Cost 

Space Shuttle 

$ 9,100 

Atlas 

$ 13,900 

Delta 

$ 15,300 


Figure 1-2 


The question becomes why is this method and research necessary and important. 
Conventional wisdom suggests that this problem has been solved. However, optimal 
control theory has demonstrated in many instances that the most intuitive solution is not 
necessarily the optimal solution. For example, to say that the shortest distance between 
two points is a straight line is not always true when one adds additional complexities or 
constraints, such as if the points are located on the surfece of the Earth, to the system. 
Monetarily, this question should be examined in detail. Figure 1-2 shows the cost of 
raising a single kilogram of material into low-Earth orbit. Any savings at all over 
conventional transfer methods would soon add up to hundreds of thousands of dollars or 
more over the life of the spacecraft. The possibility of savings alone justifies this 
examination of optimal control theory applied to this problem 
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Optimal control centers on the formulation and minimization of a cost fimction, in this 
case a function that describes the propellant expended. However, orbital motion is cyclic 
in nature. This means that the motion is periodic and repeats over a given time period. 

For many processes, periodic controls are more efficient than steady state operations. 
Acknowledging that fact, the problem then requires the examination of a further 
specialization of optimal control, namely, optimal periodic control. This process adds the 
complexity of periodic states and controls that allow a better description and 
representation of the periodic process. Through the use of numerical techniques, the 
optimal periodic control problem can be solved, creating a time history of states and 
controls that minimize the system’s cost function. 

The first step in examining this system and the orbit maintenance solution begins with 
the equations of motion, which are normalized for properly scaled numerical methods. 
Within the equations of motion, the state and control variables are identified. Boundary 
conditions are then determined with some of the conditions being a set of periodic 
functions. A spectral collocation method is used to numerically solve the optimal periodic 
control problem. This method seeks polynomial approximations for the states and 
controls in terms of their values at certain points or nodes. These nodes are called the 
Legendre-Gauss-Lobatto (LGL) points and are located at the roots of the first derivative 
of an n*-degree Legendre polynomial. The performance measure for the periodic process 
is approximated by the creation of a cost function and referenced to a steady state 
solution. Through the use of the collocation method, the OPC problem is converted to a 
nonlinear programming problem that is solved using existing routines. The full 
optimization code is then optimized for quicker performance with more reliable results by 
introducing period constraints and eliminating certain variables. An example is tracked 
throughout the thesis to highlight the different aspects of the implementation of the OPC 
and the associated code for the orbit maintenance problem. Finally, the results of the 
optimization code are compared with existing trajectories and any fuel savings are noted. 
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n. ORBITAL MOTION AND OPTIMAL CONTROL 


A. OPTIMAL CONTROL TEDEORY 

For any problem that has an infinite number of solutions, fiuther means must be used 
in order to choose the solution that maximizes the performance of the system at the lowest 
cost. The method designed to do this is called optimal control. The use of optimal 
control begins with a system that is described by a set of state variables, which describe its 
various physical limitations. Inputs into this system are called controls, which perturb the 
states within the system constraints over the time history of the problem. The 
performance of a set of controls versus another possible solution set is compared by a cost 
fimction that is minimized or maximized by the optimal control algorithm. The control 
set that minimizes the cost function is the optimal solution. Classical linear control 
systems generally relied upon "trial and error" processes that used numerous iterative 
techniques to determine the optimum set of controls. However, modem day problems 
require system dynamics that are too complex with many dififerent performance criteria to 
be met, and classical methods are insufficient to solve these types of problems. The 
complenty of the problem is increased by another magnitude when periodicity is added. 

In the optimal periodic control problem, the states, controls, and perhaps most importantly 
the boundary functions can all be periodic. Orbital motion is a periodic problem. A 
spacecraft starts from a certain location, in this case defined by its velodty, radius, and 
flight path angle. After a time period that is either given or formulated by the control 
problem, the spacecraft must return to its original states. The process begins again 
continuing in a cyclic maimer. The performance criterion used for this type of problem is 
the amount of fuel required ovct the problem's period, x. 
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B. EQUATIONS OF MOTION 


To first attack the optimal control problem the equations of motion for orbital flight 
must be obtained in order to mathematically model the physical system. The objective of 
the modeling process is to create the simplest mathematical description of the system to 
predict its response to any given input. The states are governed by a system of first order 
differential equations. 



In inertial space (N): 


^r = v 

and (1) 



m 

To find a difierential expression for the radius, state r, the derivation begins by examining 
the radius in terms of its local body coordinates. 

r = A (2) 


Taking the derivative and noting that = 0, equation (2) becomes; 

= A (3) 

To transform the expression into the inertial coordinate system the transportation theorem 
is used and equation (3) becomes: 
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w r=Bf+N(Q B X f = 5 (4) 

Simplifying: 

v = fB2+ (5) 

From Figure 2.1 another expression for v can be obtained. 

V = v-cos(r)B^ +v-sin(y)B^ (6) 

Since the expressions fort?, equations (5) and (6), are equal, the individual components 
along each body axis are also equal. An expression for f can now be determined. 

r = vsin(y) (7) 

The derivation for v follows the same general procedure as above but is worth some 
examination. The vector v is defined in its own local coordinate system, fi'ame A fi'om 

Figure 2.1. Aj is the unit vector in the direction of the velocity vector while is in the 

normal direction. Taking the derivative in the A frame; 

^5 = Ml (8) 


Using the transportation theorem to put v in terms of the inertial coordinate system, 
equation (8) becomes: 

From equation (1) and the summation of forces fi'om Figure 2.1: 

F Tcos(a)-D . . , . o . . T-sinfa) . 

m" - m - A,-g sm(r) A,-g-cos(r)-A^ +- ~-A^ ( 10 ) 

Equations (9) and (10) combine to give an expression for . 

. Tcos(a)-D 

^ =- - - i‘Sin(Y) ( 11 ) 

The above expressions for f and n are two of the five states needed to adequately 
describe orbital motion. The other states are derived in a similar manner. The full set of 
equations of motion become: 

f = v-sin(r) (7) 

. T-cos(a)-D . , ^ 

P =- - - -gstn(r) (11) 
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( 12 ) 


r = 


r 

-T 


■8 - 


cos(j') T- 8111 ( 0 ) 


m-v 


m = 


V 

0 = — -cos(y) 

r 


(13) 

(14) 


C. NORMALIZATION 

The optimal periodic control problem is not solvable by analytical means and must be 
solved numerically. A numerical nonlinear code is used to solve for the states and controls 
and determine the best solution. The size of the physical constants involved with the 
problem of orbital motion varies greatly. For example, the typical atmospheric density at 
the altitudes in question lies approximately in the lO '^ kg/m^ range while the orbital 
velocity at the same altitudes is around 10^ m/s. The numerical disparity between these 
two numbers poses significant computational difiBculty for any code. Equation (11) which 
is the first order differential of the velocity includes the drag term, D, which includes 
atmospheric density. To better control the scaling of the physical constants, non- 
dimensional units are derived through the process of normalization. 

For practical reasons, all numerical calculations are based upon reference values taken 
at a specific value in space. One of the popular ways of normalization in orbital motion is 
the creation of the canonical unit system. For Earth based systems, canonical units are 
referenced with large measurable Earth constants. One distance unit (DU) is equal to the 
radius of the Earth (6378 km). A velocity unit (VU) is defined as the orbital velocity of an 
i ma g inar y object that circles the Earth at a distance of one DU, or simply at the surfece of 
the Earth. A time unit (TU) is the amount of time it takes the same object to travel one 
radian about the orbit’s center. The gravitational parameter, yu, turns out to be one 
DU^/TU^ or simply have a numerical value of one. 

While canonical units are convenient in describing earth systems, a slight modification 
of the reference values decreases the level of numerical complexity for this problem. The 
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distance unit is referenced to the altitude where the satellite operates instead of the Earth’s 
surface. 


r = 


( 15 ) 




Numerically the code uses an altitude of300 km for the basis of all radius calculations. In 
this case, becomes approximately 6678 km. When the vehicle is at an altitude of 300 
km above the Earth’s sur&ce the value of f is one. 

Dividing the numerical value by the orbital velocity at the vehicle’s reference altitude 
normalizes the velocity parameter. 

_ V V 


v = - 


V 


ref 


•ref 


(16) 


Similarly mass and time can be normalized. 


m 


m=- 


m 


ref 


t 


f=— where = 




ref 


V 


ref 



(17) 


( 18 ) 


The first order differential equation for the radius can now be fully described with non- 
dimensional units. Equations (15) through (18) combine to allow the following derivation: 


.dr 

r = — = v-stn(r) 


Substituting the reference values: 


dr 

di 


Kef^rtf 


v-sin(Y) 


( 7 ) 


(19) 


•ref 

f -\-v-sin(y) (20) 

The structure of the first order dififerential does not change due to no rmalizati on^ rather 
the numerical values must be reinterpreted as normalized values. The normalization of the 
dififerential velocity follows the same procedure but requires an extra referencing of thrust 
and drag. 
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Breaking equation (11) into two parts: 


- 

v = — 




'ref 

where A, = 


(-g-sin(r)) 




T -co$(a)-D 


( 22 ) 


m 


The term Ai is the acceleration due to the system forces along the tangent of the orbit. 
Substituting in for the gravitational acceleration, g: 


-f. 


i> = - 


f/ 


V 


sin(x) + -A, 




=2 


V 


(23) 


'ff 


Realizing that g = —: 


v = -g.sin(r) + a, 
t. 


fl. = 


^ref 


V 


(24) 

(25) 




a. = 


tref T-cos(a)-D 


m 


Vref^ref 

Equation (26) would be sufiScient in normalizing the velocity differential "with thrust and 
drag expressed as values normalized by the centrifugal force. 

T . - D 


(26) 


T = 


/ 2'\ 

OilU JL/ — 

( *« T1 2^ 

^ref^tef 

K ^rtf ) 


[ ^ref J 


(27) 


However this form would make the relationship between thrust and drag difficult to 
visualize. Instead, thrust and drag are divided by the value of drag that the vehicle would 
experience at the reference altitude. When the vehicle is at the reference altitude the value 
of the normalized drag would be one while the value of thrust would be in terms of the 
drag force. A normalized thrust value of 10 would allow the spacecraft to have a 
thrusting force ten times the force of drag. 


(28) 


- T — D 

T = —- and D=^ 


The value for Dr^is defined by the formula for atmospheric drag. 

(29) 

The variable pr^is the density at the reference altitude, C/jis the drag Coefficient, and A is 
the cross sectional area of the vehicle along the velocity vector. Using this method of 
referencing thrust and drag equation (26) becomes: 


T ■cos(a)-D 


Kffnref 


Looking at only the coefficient of the above equation: 


Kef^ref ^ref^ref 
^ref^ref 


^refPref f C^A 

2 J \ m . 


A final variable is introduced, the ballistic coefficient of the spacecraft, B. The ballistic 
coefficient is a fimction of the vehicle’s mass and cross sectional area and is represented by 
the following equation: 


To normalize the ballistic coefficient an aihhrary reference coefficient was chosen and is 
shown below: 


^ref Pref 


Using equations (31) through (33) the following can be shown: 

Kef^ref 1 

^ref^ref ” B 

Finally, the full non-dimensional form of the velocity difterential is determined by 
combining equation (24) with the above expression. 

^ , T cos(a)-D 

v = -g- stn(r) + -==- 

mB 
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At this point all of the variables have been normalized and the remainder of the non- 
dimensionalized equations of motion are listed below. 


r = v-sin(y) 

(20) 

. . , T 005 ( 0 )-D 
v = -gsin(r)+ 

(36) 

^ fv^ '] cos(y) T-sin(a) 

(37) 

. -T 

(38) 

— V 

0 = —.cosfrj 
r 

(39) 
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ffl. ORBIT RAISING MANEUVERS 


A. THE LOW EARTH ATMOSPHERE 

While drag is considered mainly for air-breathing systems, it is the principal non- 
gravitational force that affects all objects in low-Earth-orbit. While the density of the 
atmosphere at orbital altitudes is about 10’^^ times smaller than at the surface of the Earth, 
oibital velocities are approximately 7700 m/s. Recalling equation (29) for Dref, the 
magnitude of the drag force is proportional to the square of the velocity. 



Figure 3-1 


Figure 3-1 shows the various inputs into the drag equation for the some typical low- 
Earth orbit altitudes. The variables combine to make drag a small force that continually 
acts to perturb the original orbit. By dissipating energy, drag causes the orbit to shrink. 
As the orbit altitude deceases atmospheric density increases. This increases the net drag 
force, which increases the rate of orbital energy dissipation. In the determination of 
orbital motion that includes drag, it is important to model the force of drag upon the 
vehicle at different altitudes. To acciuately describe drag, one must choose a fairly 
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accurate density model. The Earth’s thermosphere, which begins above 90 km in altitude, 
is the region of the atmosphere that absorbs extreme ultraviolet radiation. The 
temperatures of the molecules within this region vary widely between day and night due to 
different levels of UV absorption. Increases in temperature create an increase in 
atmospheric density. Other disturbances such as geomagnetic activity and solar cycle can 
vary density values on an hourly basis. Complex models have been used to describe the 
variations in density, the two most popular being the Jacchia and the Mass Spectrometer 
Incoherent Scatter (MSIS) models. Numerically, both of the models are complex and 
approximate atmospheric densities as a function of time and position. To avoid the 
complexity during implementation of the optimal periodic control problem, a simpler 
density model was chosen. The main reason for stripping away higher-order terms and 
variables in the atmospheric density model is to study the problem at the fundamental 
level. The easiest choice for an atmospheric model would be to assume a constant density 
model. When looking at oibit reboost problems the difference between upper and lower 
orbital altitudes is usually fairly small on the order of 10-50 km. The density change 
between 10 km is less than 20 percent and provides a first guess into the interaction of 
drag and optimal thrust control. While this difference may seem significant, larger 
variations occur at the transition between night and day at orbital altitudes. 

The next step in adding complexity is to assume an exponential function of altitude for 
the densities. The advantage of this model is that density closely follows the nominal 
measured results. The equation for the exponential density model is shown below. 

p = pjt h (40) 

The atmosphere is broken into separate bands with reference values given at a specific 
value taken from the MSIS model. For example, the density at 300 km is 1.87e-l 1 kg/m^ 
with the scale height, h, having a value 50.3 km. All altitudes between 275 km to 325 km 
are referenced to the values at r = 300km. A simple MATLAB function called density, 
included in Appendix 1, shows the atmospheric density determined as a function of orbital 
altitude. Figure 3.2 shows the output of the function (sohd line) with the ‘X’s 
representing the MSIS values at the given orbital altitude. The exponential model was 
used in the formulation of the optimal periodic control problem. The function density is 
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Density Model 



'-■-1-.-1_I 

150 200 250 300 350 400 

Orbital Altitude (km) 

Figure 3-2 


called during the numerical iterations to accurately define the atmospheric density at the 
local altitude. The exponential model, while increasing the complexity of the dynamical 
model, is a fairly good approximation of the density as a function of altitude. 

B. ORBIT-RAISmG MANEUVERS 

With the environment and physical laws that govern the motion of the spacecraft, the 
next step in understanding the orbit maintenance problem is to discuss orbit-raising 
techniques. Simply put, orbit maintenance is the process that maintains the satellite in a 
specific region of space during the lifetime of the satellite. As stated earlier, drag is the 
largest non-gravitational force that affects a spacecraft, specifically in terms of its orbital 
altitude. The lower altitude limit represents the lowest altitude to which an orbit can 
decay. Any further loss of altitude would signify a larger force of drag, jeopardizing the 
continued motion around the Earth. The upper altitude limit is the maximum altitude at 
which the spacecraft can operate effectively and safely. Reasons for this limit may include 
the maintaining of earth observation resolution or reduce the risk of high altitude 
radiation. Obviously of the two limits, the lower limit is the most critical due to the 
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immediate threat to the spacecraft’s lifetime. In formulating the optimal control problem 
two methods were examined to use the limits as constraints. The first type is the upper 
bound unconstrained problem. Here the satellite begins at the lower most limit and forced 
to remain above the minimum altitude. The second type includes an upper bound on the 
altitude. The difference between the two altitudes is called the orbital band. 

Figure 3-3 shows the different types of orbital transfers. The high-energy transfer is 



the tniniiniiTn time transfer but uses large amounts of propellant. The second is the most 
common method of orbit altitude maintenance and is called the Hohmann maneuver, 
designed to keep the spacecraft within the orbital band. As the vehicle approaches the 
lower altitude due to drag-induced decay, the first of two thruster bums is applied to the 
spacecraft to reboost it to the top of its band. After the first bum the vehicle travels 
approximately along an elliptical transfer orbit which has an apogee equal to the upper 
altitude limit. As it approaches apogee a second bum is made circularizing the orbit at the 
upper altitude. After the satellite’s orbit is decayed by drag the process is repeated. The 
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amount of fuel required to accomplish these bums over time is called the propeDant 
budget and directly determines the length of time a satellite will be able to maintain 
operational status. Figure 3-4 shows the typical ttohmann trajectory. The path was 
determined by a first order differential solver included in Appendix B called orbprop. 
This program takes the equations of motion with an initial condition and propagates the 
equations over a fixed period of time. During the orbital propagation an exponential 


Radius versus Orbit 



Figure 3-4 

density model is used. The actual thrusting logic sequence comes fi^om a separate 
program called hoh that is included in Appendix C. All of the Hohmaim thruster bums are 
impulsive maneuvers, meaning that the bums occur instantaneously with infinite force. 

The ideal Hohmann neglects gravity and drag loss terms. For real-world applications, a 
margin must be added to the required propellant mass to account for these losses. 
Depending on the size of the orbital reboost thmster, between five to ten percent 
additional propellant may be added to the mass budget. 
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The final two transfers of Figure 3-3 are also important to mention. First the chemical 
transfer is a real world application of the Hohmann transfer. Since impulsive bums are 
infeasible, the transfer bums are broken into distinct parts. As the satellite approaches 
perigee a finite bum is accomplished, placing the satellite into an elliptical transfer orbit. 
The apogee of this bum is less than the desired final altitude. The satellite then makes an 
additional bum at the following perigee, placing it into another elliptical orbit with a higher 
eccentricity. The process is continued until the apogee height of the transfer orbit is equal 
to the desired altitude. The total amount of bums nearly approximates the Av required by 
the Hohmann transfers. The final bum fi-om Figure 3-3 is the low thmst electrical 
propulsion transfer. Here the thmst of the electric engine is low and therefore causes a 
trajectoiy that spirals fi-om the lower altitude until the higher altitude is reached. This 
trajectory resembles the opposite of the drag induced decay spiral. 

The next maneuver to analyze is called forced Keplerian trajectories (FKT). An FKT 
is simply a thmst drag cancellation process. The thmsters are fired in opposition of the 
drag force. This requires that the thmster act in a continuous mode with the ability to 
change its thmsting force depending on the periodic variations of the local atmospheric 
density or velocity and altitude for an elliptical orbit. The trajectory would eliminate any 
altitude loss due to drag. If the orbit were circular the resultant trajectory of an FKT 
maneuver would remain circular at the ori^al altitude. Comparing this trajectory and 
propellant requirement versus a Hohmann transfer requires an additional distinction. Since 
the satellite in a Hohmann maneuver travels from the top of the orbital band to the bottom 
one must account for the variations of the drag force. Thus to cover the entire orbital 
band, three different FKT maneuvers were considered. The first, labeled the low-FKT, 
maintains the satellite at the bottom of the band. TTie mid-FKT main tains an altitude at the 
center of the band while a high-FKT keeps the satellite at the top of the band. The same 
propagation code, orbprop, tabulated the propellant used to counteract drag at the three 
different altitudes. These values were then compared to the fuel requirements of the 
Hohmann maneuver and shown in Figure 3-5. 
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0.03 


Propellant Comparison 



Figure 3-5 

The ideal Hohmaim is closely approximated to the mid-FKT. To keep a vehicle fixed at 
the lower altitude limit requires the most fuel of the four cases. The force of drag is the 
greatest while at the lower altitude. The orbital velocity is higher and the atmospheric 
density is also higher. Thus the necessary propellant to counteract the higher force must 
be hi gher than altitudes where drag is less. A common sense rule of thumb leads one to 
suggest that to decrease the amount of fuel for orbital maintenance, the satellite must be 
maintaine d at a higher altitude. Unfortunately, this reasoning may not be amenable with 
the vehicle’s mission and requirements. 

If there are no state constraints, Ross, et al [Ref 3,4], have shown that the FKT is not 
the fiiel optimal solution. This was accomplished by considering the totality of extremal 
arcs. Ignoring the special case when the maximum available thrust equals drag, the FKT is 
not a singu lar arc and thus not fuel optimal [Ref 3,4]. Since the Hohmann transfer does 
not do better than the mid-FKT the periodic Hohmann reboost cannot be the fuel optimal 
solution as well. Thus there must exist some other trajectory that is more fuel-efficient. If 
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one considers the drag loss, bum times and orbital positioning some savings of fuel can be 
made. To determine this fuel-efficient trajectory optimal periodic control theory is 
explored. 
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TV. OPTIMAL PERIODIC CONTROL 


A. PROBLEM FORMULATION 

Orbital motion in its simplest form is a cyclic process. Any orbital body subject to no 
paturbations will travel along a trajectory in a periodic manner with the orbit defined by a 
specific set of states and period. Periodicity is a further specialization of optimal control 
theory where the states and controls are cyclic and their values repeat over an optimal 
period. However before the periodic nature of the theory is examined, the general 
formulation of the optimal control problem must be established. In differential form of 
optimal control, the fisUowing represents the state and control equations. 

0 = fix(t), u(t),t) / € 9? + = [0, oo) 

(41) 

u^(t) x,(t) 

u(t)= : , x(t)= : 

yp(0\ 

Here the normalized states are represraited by the equations (20) and (36)-(39). The 
controls considered in this problem are thrust, T, and the thruster angle, a. This 
representation develops a simple mathematical form that adequately predicts the response 
of any system. The history of control inputs fi-om [to, //] is called the control history, u(-), 
while the state values over the same time interval is the system's trajectory. The trajectory 
and the control history must satisfy all of the system state and control constraints in order 
to be classified as admissible. The constraints reduce the range of values that can be 
assumed by the state and control variables. In the case of orbital motion several obvious 
constraints exist. The state equations are constrained in that the mass can never be 
negative or using a lower altitude limit places an inequality constraint upon the radius. 

The controls are constrained as well: Thrust is bounded below by zero, the force when not 
firing, and bounded above by a maximum, derived by the physical limitations of the motor. 

The next characteristic of the optimal control problem is the need for a performance 
measure. Optimal control is defined as the process which minimiz es a given performance 
criterion. In general the performance measure can be written as the following: 
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(42) 


V 

J = h(x(tf),tf)+jg(u( t), x( t),t)dt 

h 

For orbit maintenance the main performance measure is to minimize the amount of 
propellant expended. For periodic control the performance measure becomes an average 
cost measure. Under the framework of optimal periodic control theory and the specific 
problem, the performance criteria, or cost function, is simply the average amount of 
propellant required over an optimal period, z. The h portion of equation (42), can be 
eliminated through the use of judicious referencing of the state values. The periodic 
performance equation becomes: 

J^^gWMO.Odl. (43) 

To 

The functional g is dependent on functions containing the controls and the states. To 
create a fuel-optimal orbit maintenance profile, one examines the minimization of the 
propellant required over a given period. The representation of the cost function begins 
with the following equation. 

T 

Neglecting any pressure differences in the nozzle exhaust, an equation for the average 
thrust becomes: 

T=mv, = {n^ o)-m( zj)-r- v, (45) 


Putting the cost function in terms of an integral as in equation (43): 

1 f r 


^ 'T J V 


T* V 
^ 0 


(46) 


Substituting in the non-dimensional variables the periodic cost function becomes: 


(47) 


The added specialization of optimal control theory is the periodic behavior of the 
states and controls. Called optimal periodic control (OPC), the significance of this theory 
lies in its boundary conditions. In standard OPC theory, all the states are periodic. This 
results in periodic costates as well as a transversality condition, in terms of the 
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Hamiltonian, from which the optimal period may be determined. One of the goals of 
periodic control is to identify the optimal period in order to gain the most efficient system. 
However in the case of orbital motion and maintenance, the problem formulation requires 
that some of the states are not periodic. For the orbit problem, the radius, velocity, and 
flight path angle (/) are all periodic values while mass and the position angle theta are 
aperiodic. For instance mass is a state that begins at a certain value and decreases imtil it is 
zero. Angular position, 0, is also not necessarily periodic. Since some of the values are 
aperiodic while others are periodic it is difficult to apply the full optimal periodic control 
theory. Instead the problem is solved by a numaical method that is discussed further in a 
following section. An advantage of using this method is that the initial conditions along 
with the fixel-optimal trajectory are obtained. The nonlinear code is given a set of initial 
guesses for the states and controls and then is free to change the value of the variables 
assuming they are not manually fixed as a constraint. The initial conditions determined by 
the optimal code are those which minimize the cost function. 

The period, represented by x, can either be determined by the optimal code as a free 
variable or fixed. 

f(0) = r(T) f(0) = f(T) v(0)=v(t) 

m(0) = l 0(0) = 0 (48) 

m(T) = free 0(t) = free 

Equation (48) shows the boundary conditions for the case where the initial radius, velocity 
and flight path angle are all free. Variations of the boundary condition set will be 
explained with each individual case. 

A further modification was made to the cost function in order to help quantify any 
efficiency gained by the numerical code. If one assumes a low FKT trajectory, specifically 
a satellite that uses a thrust drag cancellation profile at the altitude of its lower limit, the 
cost function for this trajectory is written below. 

FKT 

T 0 v^a 

Since the trajectory is an FKT, thrust is equal to drag. Recalling equation (28), thrust is 
referenced in terms of drag at a specific reference altitude. 
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T 


D 


(28) 


T = 


D 


«/ 


and 



Using this fact, an FKT trajectory would have a value T = 1. Equation (49) becomes: 

Using equation (50) as a reference one can then divide the cost function obtained in 
equation (47) to have a referenced performance measure. 

Jv 1 

J = -^ = -]Tdt (51) 

J Fkt ^ ® 

This referenced cost fiinction is very relevant when the radius is constrained to never drop 
below the lower altitude limit. A cost function with a value of one would equal the 
performance of the low-FKT. A cost function less then one would mean that the optimal 
trajectory would be more efficient than a low-FKT. In other words, the aim of the code is 
to determine a J less than one. 


B. NUMERICAL METHOD 

Traditionally, optimal control problems are solved using shooting methods that require 
the formulation of costate equations obtained fi'om first order necessary optimality 
concfitions. There exist no widely accepted methods and extensive theoretical results for 
determining the necessary conditions for a partially periodic problem. Therefore the 
emphasis of this thesis is to solve the optimal periodic control problem directly. A speptral 
collocation method is used to numerically solve the OPC problem for the proposed states 
and controls. The spectral collocation method used here has already been successfully 
applied to solve a class of linear and nonlinear optimal control problems >vith state an^ 
control constraints [Ref 5,6]. 

The premise of this method is to discretize the nonlinear control problem at specific 
points or nodes and convert it into a system of nonlinear algebraic equations with 
unknovm as the values of the states and controls at the nodes. The resulting nonlinear 
program can be solved by existing routines. In order to create orthogonal polynomial 
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approximations for the control and state equations a spectral collation method is e mp lnyp-d 
to determine the values of these functions at the Legendre-Gauss-Lobatto (LGL) points. 
These points are used as the collocation points and Lagrange polynomials are used as 
orthogonal trial functions. 

Given any function F(t) that exists over an interval from [0 ,t] a time transformation 
must be used to convert it into the interval [-1,1], where the Legendre-Gauss-Lobatto 
points lie. 

where (52) 

te[0,r],- /e[-l, l] 

Using (43) and (52) the performance function becomes: 

1 1 

J^-\g(x(t),u(t),t)di (53) 

And the periodic boundary conditions change to: 

x(-l) = x(l) (54) 

Let Ln( i ) represent the Legendre polynomial of order N. These polynomials are 
determined by the following recmrsive expression: 

(i + (t) - (2i + \)tL, (t) + (t) = 0 

where (t) = 1, (f; = t, (55) 


and i = \X---,N 

Existing numerical codes were used to determine the location of the zeros that give the 
time interval between the nodes. The t i, where i = 0,1.. .N, are defined asfo = -l,tN = 
1, and t k equal to the zeros of L^(t ). The next step is to construct the polynomial 
approximations by first defining the Lagrange polynomials in terms of Ln( t )s. 




1 ((=-i)4('0 

+ i-t, 


(56) 


It can be shown that: 
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Using the above relations, one can now build the polynomial approximation of the 
following form: 

y=0 


/ = u...,p 

y=o 


By differentiating equation (58) and introducing the differentiation matrix Djk we get the 
following ejq)ression: 

i,''(ij)=ip^x,(t,) (61 


The differentiation matrix Djk is defined as follows: 

' LJh) 1 
-n{n+\) 


Djk- 


4 

Ar(Ar+i) 

4 

0 


jz=k = 0 

j = k = N 
otherwise 


Placing the equations (59) and (60) in vector form: 


Jt=0 

*=0 

i=0 

where 

Ok ~ [^\k mk ] 

bk — [b^k’^2k>"-^pk] 


The coefficients a* and bk are yet to be determined but it should be noted that: 

b^=u^(tk) 
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The next step in the approximation process is to discretize the cost function integral. It 
can be shown that the cost function can be rewritten in the following form. 

'Zg(a„b„ 4 )w, (64) 

Z it=0 


Where the term Wk are the weights given at each node and are defined as; 

2 1 


(65) 


The state equations are discretized by substitution and collocation at the LGL nodes, . 

The equations are then converted into the following series of algebraic equations in terms 
of the vectors ay and bj. 




N 

where d,^ = 

y=o 


( 66 ) 


and k = 0,...,N 

Similarly the system constraints can be approximated in the same manner. 

B,=g{x^(i,),u^(t,))<0, 


(67) 


C. NONLINEAR PROGRAMMING CODE 

The next section will describe the various computer codes that were built in order to 
solve the optimal control problem directly. Each program represents a single or series of 
functions. When used in conjunction with each other the result is a series of commands 
which solves the optimal periodic control problem, checks the values against a differential 
solver, compares the results against different trajectories and finally records the data. All 
the codes are programmed m MATLAB and it is assumed that the reader has a 
backgroimd in the MATLAB language 
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1. Orhprop - Appendix B 


Mentioned in the previous section, orbprop is a first order differential equation 
solver which propagates the equations of motion given a set of initial conditions. Figure 
4.1 shows the declarations for the five states. 


% States 

% x(l) = radius 

% x(2) = velocity 

% x(3) = angle gamma 

% x(4) = mass 

% x(5) = angle theta 


Figure 4-1 

Orbprop is a function that is called by other functions. The program requires that initial 
conditions be given for each of the states and time. The program then takes a profile for 
the controls, thrust and thruster angle, and chooses time steps in order to begin 
propagating the values. The profile for the controls depends on which program calls 
orbprop into action. For instance, the results fi'om the optimal code can be transferred 
into a time control history that this program can use to propagate the states. Other cases 
include setting the thrust, denoted by Tp in Figure 4-2, to equal drag for the FKT 
trajectory while another trajectory would call for thrust to equal zero allowing orbital 
decay in a “free faD” manner. Figure 4-2 is an excerpt fi'om the program that shows the 
diflferential equations of motion. Tp and alphap represent the interpolated values of the 
thrust and thruster angle respectively. Since the propagator uses different time steps than 
the location of the LGL points, interpolation is required to define the controls at these new 
points. 
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As = (Tp*cos(alphap)-D)/(x(4)*B); 

An = Tp*sin(dphapy(x(4)*B); 

ddot(l) = x(2)*sin(x(3)); 

ddot(2) = -g*sin(x(3))+As; 

ddot(3) = (x(2)^2/x(l)-g)*(cos(x(3))/x(2))+An/x(2); 

ddot(4) = -abs(Tp)/(ve*B); 

ddot(5) = x(2)*(cos(x(3))/x(l)); 

Figure 4-2 

The expressions for As and An in the figure represent the only terms in which the controls 
affect the equations of motion. The five “ddot” expressions are the normalized equations 
of motion derived earlier (equations (20),(36)-(39)). 

2. Appendix D 

This program is the macro program that controls the flow of the optimal periodic 
control problem. It begins with a series of variable declarations that give the problem 
physical meaning and definition. The next step calculates the normalized ballistic 
coefficient. This term is the only term that relates to the actual spacecraft being modeled. 
Its mass, area, and coeffident of drag are all contained within this single term. A very 
large value for the ballistic coefficient means that the spacecraft is either very massive and 
thus not greatly affected by drag or that the spacecraft is small and drag resistant. Low 
values of B would suggest that the spacecraft is very susceptible to drag due to its small 
mass and/or large surface area. 

The next section of this program initializes the vector that contains all of the states and 
controls, named “aop”. Figure 4-3 calls a MATLAB program included in the 
Optimization Toolbox called constr. This function finds the constrained minimum of a 
function of several variables. 

aop=constr('orbcrit',aop,options,vlb,vub); 

Figure 4-3 
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This command calls the function file called orbcrit, which is discussed in the following 
section. The options are a set of choices the user can alter to customize the minimization 
process including setting the minimization tolerances, number of iterations, etc. The “vlb” 
and “vub” are the optional lower and upper bounds. In most cases these were not used, 
instead all altitude constraints were entered in as inequality constrdnt equations. After the 


i=aop(l:n); 
v = aop(n+l:2*n); 
gamma = aop(2*n+l:3*n); 
mass = aop(3*n+l:4*n); 

T = aop(4*n+l:5*n); 
alpha = aop(5*n+l:6*n); 
theta = aop(6*n+l:7*n); 


Figure 4-4 

successful completion of the optimization routine the vector “aop” was defined into all of 
the states and controls shown in Figure 4-4. During this final declaration the variables 
represent the optimal states and controls which minimize the given cost function. These 
variables will later be compared to different trajectories in order to measure the 
effectiveness of the minimization code. 


3. Orbcrit - Appendix E 

This program is the center of the spectral collocation method. The code begins 
with the creation of the cost function. Figure 4-5 contains a section of orbcrit that 

for i=l:n 
fii(i)=aop(4*n+i); 
end; 

costfii= l/(2)*sum(w.*fh'); 

Figure 4-5 

pertains to the creation of the cost function. The vector “fii” is given the control history 
values of the thrust and is then summed in the “costfh” expression. Recalling equation 
(51), the derived expression for the reference cost function; 
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The cost function is simply the sum of the weighted values of the thrust. 

The next part in the orbcrit code is the creation of the state constraint equations. 
Figure 4.6 shows the five state constraint equations. 


% Radius 

g(i)= (2/tau)*sum(Dn(i,;). *aop(l :n))-aop(n+i)*sm(aop(2*n+i)); 

% Velocity 

g(i+n)= (2/tau)*sum(Dn(i,:). *aop(n+l :2*n))+(G(i)*sin(aop(2*n+i)))-Ast; 

% Gamma 

g(i+2*n)= (2/tau)*sum(Dn(i,:).*aop(2*n+l :3*n)); 
g(i+2*n^g(i+2*n)-(((aop(n+i)^2)/aop(i))-G(i))*(cos(aop(2*n-i-i))/aop(i+n)); 
g(i+2*n^ g(i+2*n)- Ant/aop(n+i); 

%Mass 

g(i+3 *n)= (2/tau)*sum(Dn(i,:). *aop(3 *n+l :4*n))+(aop(4*n+i)/(ve*B)); 

% Theta 

g(i+4*n) = (2/tau)*sum(Dn(i,:).*aop(6*n+l :7*n))-aop(n+i)/aop(i)*cos(aop(2*n-+i))); 


Figure 4-6 


Before discussing Figure 4-6 in detail a brief explanation is required concerning the 
way in which constraints are implemented for the optimization routine. All constraints are 
referenced to zero. For example, to place a constraint on the initial radius the first element 
of the radius vector is modified. 

r(tj = \ (68) 

Recalling that the computer code uses the vector “aop” to represent the states and 
controls and that radius occupies the first n terms (where n equals the number of LGL 
points) of the “aop” matrix, equation (68) can be referenced to zero and put in the correct 
form. 

aop(\) -1=0 (69) 

The process for inequality constraints is similar. In order to place an upper bound on the 
amoimt of thrust a given engine can produce, the following expression would be used. 

T<5 


(70) 





Figure 4-5 shows the part of the “aop” matrix that relates to thrust. Equation (70) 
becomes: 

aop(A* n-\-i)-5<Q where/ = L..w (71) 

Equation (71) prevents any terms within the time history of thrust to be above 5 
normalized units of thrust. 

Figure 4-6 uses the format required by MATLAB and creates 5n equality constraints 
for the five states. Taking the first equation, the radius state equation, allows closer 
examination of the constraint construction. Recalling equations (66) and (20): 

where d, = (66) 

j=0 

and k = 0,...,N 

f = v • sin(y) (20) 

Substituting equation (20) defined at the individual time steps, 4, into equation (66) the 
constraint equation becomes: 

Duff ti)-y(tk) sin( r(t^)) = 0 (72) 

T 1=0 

From Figure 4-4 the expressions for the states can be substituted in equation (72) to 
transform the equation in terms of the “aop” vector. 

2 AT 

- ^opfn + AJ ■ sin('aop('2n + kJJ = 0 (73) 

T 1=0 


The first term, (2/tau)*sum(Dn(i,;)), follows the spectral collection method where “Dn” is 
the differentiation matrix governed by equation (61) shown below. 

f 1 




^N(^k) ^k 

-NiN + l) 


4 

N{N + \) 
4 

0 


j^k 


j = k = 0 


(61) 


j = k = N 
otherwise 
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The values of the radius are continually iterated as the program is called repeatedly by the 
constr MATLAB function. The other four state constraint equations all are the same 
format, each containing the equation of motion for the particular state. 

The next two sections of orbcrit define the physical constraints and periodicity of the 
orbit maintenance problem. The remaining equality constraints are shown in Figure 4-7. 


Periodic Constraints 

g(5*n+l) = aop(l)-aop(n); r 

g(5*n+4) = aop(n+l)-aop(2*n); v 

g(5*n+2) = aop(2*n+l)-aop(3*n); y 

Aperiodic Constraints 

g(5*n+5) = aop(3*n+l)-l; m 

g(5*n+3) = aop(6*n+l); 0 


Figure 4-7 


Only the final two expressions fi^om Figure 4-7 are given fixed values. One of the equality 
constraints pertains to the initial value of the angle theta, which is set to zero. Mass is 
initially set to one through the last equality constraint of Figure 4-7. The other three 
equality expressions allow radius, velocity, and gamma to be periodic variables; the initial 
value of the state is equal to its final value. First is a series of inequality constr aint s shown 
in Figure 4-8. 


g(i+5*n+5) = - aop(4*n+i); 
g(6*n+5+i) = aop(4*n+i)-5; 

Figure 4-8 

In the most basic formulation of the problem, the only inequality constraint consid^ed 
was the upper and lower bound of the thrust. While most modem day thmsters are either 
full on or full ofl^ thrust was modeled as an engine capable of a full range of throttling 
values, ranging fi-om oflF, a value of 0, to full on, in this case a value of 5. 

Soon after the code was initially tested a new constraint was added to the Kstm Figure 
4-8. This constraint forced the initial value (and hence its final valued crforbkal radius to 
be a predefined location in space, normally at thq reference altitude. Thei reason fw this 
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constraint can best be explained as follows. The cost function, as noted before, relies on 
the amount of thrust used over the optimal period. Intuitively, the amount of drag at 
higher altitudes is less than the force of drag at lower altitudes. In an attempt to lower the 
cost function the code tries to force the radius, initial and final to be as high as possible. 
Given enough computing power and an infinite number of minimization iterations the 
radius should be infinity where the atmospheric density is zero, allowing the spacecraft to 
orbit the Earth indefinitely at no cost. Thus to study the effectiveness of this code with a 
physically feasible problem, the radius was set to an initial value. Consequently, since the 
radius is a periodic variable the final value was also set to equal the initial altitude. 

4. Other Programs - Appendices F and Beyond 

A series of additional programs used to study aspects of the fuel-optimal problem 
are included in the appendices. Several of these programs deal directly with the spectral 
collocation method. The code diffm creates the differentiation matrix required by the 
coUocation method and that is used in the program orbcrit. The programs lobatto, mxt, 
mxtj, and tqr, written by Professor Bill Gragg of the Naval Postgraduate School, 
computes the abscissa and weights for the n-point LGL quadrature problem. 

Variations of the orbcrit and orbopt programs are included. During the course of this 
thesis, modifications to the constraint set or initial state and control guess created a series 
of modified programs. The different programs will be discussed in the following sections, 
which discuss several aspects of the fuel-optimal trajectories. 

The final group of programs is the programs that use the orbprop code and are used to 
compare the trajectories derived by the optimization process versus steady state 
propagated trajectories. 
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V. DISCUSSION OF RESULTS 


A. INTRODUCTION 

While attempting to find a fuel-optimal orbit maintenance trajectory was the main goal 
of this research, a secondary goal was to test the implementation of the spectral 
collocation method through the use of non-linear programming to directly solve the 
optimal periodic control (OPC) problem. The first step was to build a fi-amework that 
contained the computer codes required to solve this type of problem in accordance with 
the OPC theory. The most basic codes were discussed in the previous chapter. To test 
the series of programs, the most general case of the fuel-optimal problem was attempted. 
The results were used to identify several trends. Using these conclusions, the general 
programs were modified by either changing the system constraints, boundary conditions, 
or reducing the number of fi^ee variables. As a result of these modifications, the orbit 
maintenance problem was studied more realistically by defining physical limits and 
constraints of current day systems. For example, the codes were modified to include an 
upper altitude limit. This constraint causes the trajectory to be restricted within an orbital 
band. Satellites that possess limited communications ability or those th^ depended on 
imagery resolution are examples of an altitude-limited platform. Each of the spedal cases 
examined are discussed in detail in following sections. 


B. TBQE FREE CASE 


The most general problem formulation allows defining the satellites physical 
characteristics, an initial altitude, and a set of initial guesses for the state and control 
histories. The bulk of the satellite’s characteristics is included in the equation for the 
normalized ballistic coefficient. 


' rrefPref^ 
2 > 


(33) 


Recalling equation (32); 
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Using equations (32) and (33) the value of the normalized balhstic coefficient is shown 
below; 


B = 


2 m 


(74) 


The reference altitude and density are predefined values and are mutually dependent. The 
spacecraft physical properties are included in the following ratio. 


Physical Characteristics = (75) 

In order to highlight the effects of drag a satellite with a low ratio of mass to area is 
desired. The reason behind this idea is that a difference in propellant usage under varying 
drag conditions is a small quantity. By increasing the losses due to drag the value of one 
trajectory can be more easily identified against another. For purposes of this work, a 
generic satellite was created that leveraged the fact that drag was the major non¬ 
conservative perturbing force. 


Generic Spacecraft 

Physical Characteristics 

Spacecraft Area 

500 m^ 

Spacecraft Mass 

3000 kg 

Maximum Thrust 

3.5 N(f =5.0) 

Altitude 

300 km 

Density 

1.87*10'" kg/m^ 

Coefficient of Drag 

2.35 


Figure 5-1 


These values are for a theoretical spacecraft but are reasonable for certain kind of 
platforms, such as space-based radars, large antenna platforms, or inflatables. Figure 5-2 
















shows the normalized ballistic coefiBcient for several different systems assuming a 
reference altitude of300 km. 


Ballistic Coefficient 

Generic Spacecraft 

4.09*10'* 

ISS-DACT6 

1.26*10® 

Space Telescope 

4.72*10® 

Landsat-1 

4.04*10® 

Echo-1 (comms) 

8.24*10^ 


Figure 5-2 


The other characteristic of the generic spacecraft is a measure of the effectiveness of 
its control, or simply the size of its orbital maintenance thruster. Defined in the 
normalization section, thrust is placed in terms of the force of drag at the reference 
altitude. Recalling equation (29): 

(29) 

Udng the values for the generic spacecraft, would equal approximately 0.7 N. 

Current orbit maintenance thrusters have thrusts ranging from imder 1 N to over 400 N.^ 
In normalized terms a thrust equal to one would be 0.7 N. The thruster modeled in the 
generic spacecraft has a maximum of 5 or about 3.5 N available. Other thruster sizes were 
examined and the results will be discussed in following sections. 

For the free case, the thruster was allowed to fire in any direction relative to the 
instantaneous satellite velocity vector. The only other constraint placed upon the system 
was that the initial altitude was fixed to the reference altitude of300 km. Since the 
spectral collocation method requires the user to input a series of initial guesses, a series of 
initial values were placed into the “aop” vector. Unfortunately, initial guesses do have 
some influence on the results obtained fix>m the optimal control problem. Thus after 
judicious ejq)erimentation the initial guesses were made as given in figure 5-3. 
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Initial guesses for the Optimization Code 


Physical Value 

Normalized Value 

Radius 

6678 km 

1 

Velocity 

7769 m/s 

1 

Mass 

3000 kg 

1 

Gamma 

1.1 radians 

1.1 

Thrust 

3.5 N 

5.0 

Tau 

720 minutes 

50 TU 

All Others 

0 

0 


Figure 5-v-3 


W^th all the information entered into the code, the optimal control program, orbopt is 
run within MATLAB. The code is run primarily on an Intel Pentium processor operating 
at 233 MHz. Total run time for this case was approximately 2 hours. The constr.m code 
in MATLAB takes the initial guesses and calculates the given cost function. It then 
creates a set of gradient information and chooses new values for the controls and states 
and recalculates the cost function. If the cost function is lower than the previous value, 
the next iteration occurs. The process continues until all of the state and control 
tolerances are met with the minimum cost function. 

The results are tabulated by a program called orbres, included in the appendix, and 
then plotted by the program orbconv. Figure 5-4 shows the resulting plots for 4 of the 5 
states versus time, which is proportional to the fifth state the angle theta. 

Figure 5-4 shows several points worthy of notice. The spacecraft is boosted to a 
higher altitude and then allowed to decay to the original altitude. The flight path angle. 
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Figure 5-5 
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the orbit is nearly circular throughout the transfer. Mass is expended rapidly during the 
first few orbits and then is constant throughout the coast/decay portion. Figure 5-5 shows 
a plot of the two controls of the spacecraft, the thrust and the thruster angle, a. 

The majority of all thrust is accomplished in the first 3 orbits of the spacecraft 
(approximately 20 radians) with an additional bum at the end of the period. The 
important point to notice is that the thmst profile is not a finite-Hohmann transfer bum. 
The optimal thmst curve dips and peaks meaning that it is not a bang-bang solution. 

In addition, unlike a Hohmann transfer, the optimal trajectory allows the thruster angle, 
measured in the variable alpha, to change. The value of alpha range form -0.1436 to 
0.3854 radians while the thruster is firing. Figure 5-6 shows graphically the variation of 
alpha during the transfer bum. 

T his range of thruster cant angles is a normal departure from cimrent orbital transfer 

Variation of Alpha 



maneuvers, which xise bums in the direction of the velocity vector only. However when 
the thruster is off, (i.e. T = 0), the thruster angle wanders between its upper and lower 
limits. This example used ± 180° as the bounds for this value. The angle at which the 
thruster points is only cmcial during engine firing. Figure 5-5 shows the wandering of 
alpha during the decay portion of the simulation. Alpha returns to the neighborhood of 
zero during the small thruster firing at the end of the run. 

The final bum is also an important characteristic of the trajectory. Recall that three of 
the states; radius, velocity, and gamma, are periodic. The initial values must equal the 
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final values. As the satellite decays, it approaches the initial altitude and velocity. 

However its angle ganama, measured between the velocity vector to a line tangential to the 
radius vector, is negative or tending to point slightly towards the Earth. The final bum is a 
correction to the angle gamma. The bum rotates the velocity vector to match its original 
direction. 

The next result of this simulation is to examine the cost function. In actuality, the cost 
function is a ratio between the fuel required for the optimal trajectory and the fuel required 
for an FKT profile at the starting altitude. Equation (51) shows the ratio mathematically. 

Jr, 1 ^- 

J = -f - \Tdt (51) 

Any cost function that is less than one means that the trajectory used less fuel than the 
FKT trajectory. For the case just considered, the cost function equals 0.784. This means 
that the optimal trajectory uses 22% less fuel than a trajectory that would follow a thmst- 
drag cancellation path at the original altitude. Taking a closer look at the physical 
characteristics of the trajectories helps explain the large difference in fuel. A satellite that 
is at a higher altitude than another will ejq)erience lower drag. The lower the drag, the 
less the fuel expended to keep the vehicle flying. Density and orbital velocity decrease 
with increasing orbital altitude. Acknowledging this physical characteristic, a first guess 
on trajectory design would be to boost the satellite as high as possible and let it decay to 
its original altitude. The optimal trajectory resembles that concept. Figure 5-3 shows the 
radius increasing to a maximxun after 3 orbits and then decay towards the original altitude. 
Thus with this is mind one can examine the trajectory more closely and determine during 
vdiich parts of the flight path might hold savings over traditional orbit maneuvers. The 
main difference lies in the initial bum and wiH be explored further in following sections. 

The next parameter to examine fi'om the simulation results is the value for the optimal 
period, x. The value of tau for this example is approximately 110 normalized time units, 
or just over one day. In a perfect simulation, one would expect that the optimal period 
would mean that the satellite must follow the trajectory repeatedly for the lowest fuel orbit 
maintenance consumptioa This idea that the period is a definable value is taken from 
optimal periodic control theory and its application to air-breathing platforms. For 
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atmospheric flight, periodic controls generate trajectories that are lower in cost than 
steady state continuous control. It is a reasonable assumption that tWs would also hold for 
atmospheric space flight. However there exists a major difference between the two 
regimes. Thruster performance in low Earth orbit is nearly independent of orbital altitude 
while air-breathing engjne performance relies heavily on atmospheric density, oxygen 
content and temperature. Drag forces are lowered with higher altitudes in the atmosphere 
but as altitude increases, engine performance decreases. At orbital altitudes the benefit of 
higher altitudes and lower drag is not offset by a decrease in engine performance. Thus 
there is no optimal point in space for the minimum trajectory. This suggests that the 
optimum altitude to fly would numerically be at infinity. 

Examining the optimal tau from this example one asks why is there an optimal period 
at all. If the optimal orbital altitude were infinity in order to achieve the fuel-optimal 
trajectory the satellite would need to be boosted to an infinite altitude and allowed to 
decay for an infinite time. In actuality, the amount that the satellite can be boosted is 
limited by the amoimt of propellant mass the satellite carries. The generic spacecraft is a 
3000 kg vehicle with perhaps 40 percent propellant mass. With 1200 kg of propellant the 
spacecraft would optimally be boosted as high as possible. The resulting period for the 
single bum and subsequent decay would be on the order of years up to infinity. 
Unfortunately the computer with its associated non-linear programming code has a 
difficult time modeling infinity. The program runs in an iterative method beginning with 
the initial guesses and recalculating the cost function after the gradients have been 
computed. The initial guess for the generic satellite example was 50 time units. The code 
increased tau to 112 time units. The code increases the period in order to lower the cost 
function since a longer period allows for a higher boost and subsequent decay. 
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However to prove that this value is the optimal period for this spacecraft and altitude, 
the code should be able to reproduce the same states and controls given a different set of 
initial guesses. For a second run the code was given the same initial conditions as the 
original example excluding the value of tau which was set at a value of 112 time units, the 
original example’s optimal tau. The simulation took just over 2 hours to run and gave a 
new optimal tau of 5513 time units, about 56 days. Thus one can conclude that the 
original example did not give the optimal tau. 

The cost function for this second run is 0.1801, or in terms of propellant, 18% of the 



fuel required for an FKT trajectory at the original altitude over the determined period. 
Again the code boosts the spacecraft to a high altitude and then allows a decay to the 
original height. In the second run the maximum altitude was approximately 1.023 distance 
units (6831 km) or about 154 km above the original altitude. The first run maximum 
altitude change was approximately 23 km. The larger altitude change allows a lower cost 
function over a longer period. Figure 5-7 shows the altitude profile for the second run. 
One of the problems with the second run is the number of LGL points conpared with the 
length of the simulation. For both runs the minimization program used 24 LGL points. At 
each of the points the states and controls are known but in the second example only 24 
known points are spread out over 5500 tune units. The simple correction for this would 
be to increase the number of nodes. However this method has difficulties with high values 
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of “n”. The higher complexity and larger matrices have the opposite desired effect. 
Computing inaccuracies and matrix scaling problems occur. These problems will be 
discussed more folly in following sections. The lack of detailed representation in the states 
and controls is one of the reasons the two graphs for the radius versus theta are different 
between the two runs. Even with the differences, both graphs imply that the trajectory 
with the minimum cost requires an initial boost to a higher altitude followed by a coast 
period towards the original altitude. Figure 5-8 shows the thrust versus theta for the 



Figure 5~8 


second run. The two graphs for the thrust are more similar. The initial thrusts are large 
with two discemable peaks. While the second run has more thruster activity after the 
initial boost, the thrust is usually small when compared to the initial bums. 

Numerically the program iterates upon the initial guess and continues imfo the 
constraints on the state and control variables are within given tolerances and when the set 
of perturbing gradients are below a given value. When guesses are sufficiently away from 
the optimum values these tolerances may be satisfied at points of local minimums. Figure 
5-9 is an attempt to visualize this conclusion and does not represent the actual plot of the 
system’s space. 
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MANY LOCAL MINIMUMS 



Figure 5-9 


Numerically the code approaches the local minimum close to the initial guess. If the 
gradients are small enough the cost hmction on dther side of the minimum is higher. If 
the ititnimum is a local minimum then the code focuses on that location ignoring the 
possibility of a global minimum location. 

If the procedure were to be repeated with the substitution of 5513 time units entered 
for the initial period the result would be an optimal period of several magnitudes higher. 
Given enough computing time and power the code would approach infinity as new states 
and controls are determined to produce a trajectory that would send a spacecraft as high 
as possible in a ringle bum to escape the effects of atmospheric drag. The next step in the 
search for the fuel-optimal trajectories is to fix the orbit maintenance period. When the 
value of tau is fi’ee the code tends to look towards infinity. The reasonable approach is 
then to look at fixed periods of time to determine if the manner in which bums are 
accomplished can lead to a one solution that is better than anoth^. This type of analyris 
can be extended to many different missions that are used today. For example, the 
Intemational Space Station (ISS) has an orbit maintenance plan that calls for the station to 
be at a specific altitude in order to rendezvous with the Space Shuttle at specific dates. 
The ISS planners calculate how high the station is to be reboosted dependant on the time 
between rendezvous. The next section vrill examine the optimal trajectories for a given 
fixed orbit maintenance period. 
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c. 


FIXED TAU 


This section takes a look at the results of the optimization code when a fixed value for 
the orbital maintenance period is used. In actuality, this section will deal with the 
numerical accuracy and effectiveness of minimization by the non-linear programming 
code. Recall that the second run of the previous section had a very low cost fimction but 
the plot for the radius did not have a smooth first bum but allowed numerous interior 
bums and decays. The following table shows the masses and periods of the first two runs. 



Check 

Program 

First Run 

Second Run 

Period 

3014.3 

112.6 

5513.4 

Propellant 

Used 

0.0311 

0.0057 

0.0637 


Figure 5-10 


The column labeled "'■Check program” is another code that estimates the fuel required for a 
continuous bum to the desired altitude followed by the time to decay to the original 
altitude. The numbers shown in this column in figure 5-10 reflect a bum to the same 
maximum altitude of the second run. This program will be further explained in the 
following sections where the optimality of the trajectories is discussed. Using the values 
in figure 5-10, a propellant consumption versus time plot is included in figure 5-11. The 
periods for the three different trajectory calculations are different; therefore, by extending 
the overall comparison period by extrapolation, a clearer picture of which trajectory is the 
most efficient is obtained. 
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Propellant Consumption 



Figure 5-11 


To determine the difference over time each of the different trajectories were repeated. For 
example when the trajectory of the &st run reached a radius equal to the original altitude 
the cycle was repeated. The difference in propellant consumption between the first and 
second runs shows the value of boosting to a higher altitude. Fuel consumption is 
dramatically reduced by a large reboost assuming the spacecraft does not have to be at a 
given altitude at frequent intervals. The small diSerence between the second run and the 
check program is harder to explain. The check program completes a constant bum at 
maximum thmst up to an altitude equal to the maximum of the second run, a radius of 
1.023 distance units. Mentioned earlier, the second run is not a smooth bum and uses 
bums during the time normally associated with the decay period. A function that would 
have contained a single bum to a higher altitude would even result in a lower cost 
function. For example if the spacecraft would have been boosted an additional 15 km 
(1.025 distance units) the mass used would have been 0.0336 units with a period of 
4261.6 TU. This period approaches the period of the second run with lower propellant 
consumption. Figure 5-12 shows the comparison with a higher altitude boost. 

The high boost (radius = 1.025 DU) is clearly the trajectory that would possess the lowest 
cost function. The question becomes as to why did the optimum code not give an 
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optimum trajectory. The answer lies in the complexity and size of the vectors involved 



Figure 5-12 


with the spectral collocation method, namely, the creation of the state and control matrices 
and the number of system constraints. 

A method to test the optimality of the non-linear code is to run the simulation for a 
fixed period equal to the period of the firee case and compare the results and cost 
fimctions. The first run of the jfree section will be used as the benchmark against which to 
compare the fixed period simulations. The first nm used the fiill equations of motion and 
all of the system constraints and gave a very smooth trajectory with a low cost function. 
The optimal period was 112.6 TU with a cost flmction of 0.7837. It will soon be obvious 
that this run is a very unique run in that the combination of guesses for the states, controls, 
and initial period all allowed the code to iterate to a near optimal solution. In 
experimentation with the code these “good” runs were very infi'equent vdth problems 
ranging fi'om matrix scaling problems to solutions that would seem to settle in a local 
minimum vice the global minimum. In any type of numerical study, reproducibility of 
results is a key tenet of success. The following runs show the ways the original optimal 
code was modified in order to reproduce the results of the original run. 
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The first step in this process is to fix the period at the optimal value obtained form the 
first run. The codes orhopt and orbcrit were modified and renamed tfopt and tfcrit, 
included in the appendix. The main difference in the programs was the removal of x fi'om 
the “aop” matrix. 

aop(7*n+l) =150; from orboptm 
tau = aop(7*n+l); 

becomes 

tau = 112.6 

Figure 5-13 

The change to the orbcrit program is similar. Figure 5-14 shows the four states as a 
result fi-om the period fixed case. The plot of the radius in the top left comer of Figure 5- 



X lO"^ Gamma vs Time Mass vs Time 



Figure 5-14 

14 is very important. Initially the code allows the radius to fall below the initial altitude of 
300 km. The general rule of thumb established in previous sections is that fuel 
consumption will decrease with high altitudes. Conversely, one would e 3 q)ect that lower 
altitudes require more propellant given the same period. This is exactly the case for this 
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nin. The cost function for the fixed period case is 0.8559 or 9.2% worse than the original 
run. Figure 5-15 shows the controls for this trajectory. 

Thmst vs Time 



Figure 5~15 

The main diflFerence between the control plots in the two nms is the position of the main 
thruster activity. In the ori^al run the majority of the thrusting action occurred within 
the first 20 time units. Figure 5-15 shows that the main reboost thrusting is not begun 
until approximately 20 TU. The initial thrust spike at time 10 is the reboost maneuver 
that takes the spacecraft that has descended below the origina! altitude back to its initial 
height. The final bum is approximately equal to the original run and is required to match 
the initial and final conditions for the periodic states. The controls for the alpha are fturly 
analogous to the original run’s thmst angles. When thmsters are firing the thmster angles 
are in the neighborhood of zero, meaning that thrust is usually directed towards the 
velocity vector. 

The results fi’om the fixed period run were disappointing in terms of result 
reprodudbility. The next step was to reduce the number of variables in the non-linear 
code further to reduce the numerical complexity of the system. The next simplification 
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affected the angle theta within the equations of motion. Recalling equations (20) and 
(36)-(39): 


r = v-sin(y) 


(20) 

^ . T ‘COs(a)-D 

v=-g-stn(r)+ - 

mB 

(36) 


cos(y) T 'Sin(a) 

V ni’V-B 

(37) 

. -f 

m- ^ 


(38) 

— V 

6 = 

r 


(39) 


The equation for 0 , equation (39), determines the differential change of theta over time. 
This value is not used in any of the other equation of motions and is conq>letely 
independent of all of the other states. Thus the removal of theta form the system of 
equations would not eflfect the other states. The loss of theta from the results does pose 
some minor problems in analyzing the data but since die values of gamma are very small, 
on the order of lO"^ radians, the angle theta can be approximated by the simulation time. 
The error in this approximation is significantly less than 1% with an orbit maintenance 
period of approximately 100. This error only changes the angular position of the 
spacecraft at any given time and does not effect the magnitude of the orbital radius. The 
codes tfopt and tfcrit were modified as shown in figure 5-16 and renamed notopt and 


Changes to orbopt. 

Removed - aop(6*n+l :7*n) = zeros(n,l); 

Reduction in the number of constraints - 
options(13)=5*n+6 to =4*n+5; 

Changes to orbcrit: 

Removed the theta equation of motion - 
g(i44*n) = (2/tau)*sum(Dn(i,:).* 
aop(6*n+l :7*n))-(aop(n+Q/aop(i)* 
cos(aop(2*n+i))); 

Removed constraint - g(5*n+31 = aoDr6*n+l): 

Figure 5-16 
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notcrit, both included in the appendix. The effect of these changes reduces the size of the 
“aop” vector by removing “n” quantities. The removal of constraints also reduces the size 
of the matrices used by MATLAB. The results reflect the changes with a more optimal 
run. Figure 5-17 shows the states of the no theta run. 


Radius vs Time 


Velocity vs Time 






Figure 5-17 

Figure 5-17 shows plots that are very similar to the original run. The radius is only 
significantly different in the first few time units. Here the radius is allowed to follow and 
descend slightly below the original altitude for a very short period. The plot for gamma is 
consistent with earlier results with a slight increase in the angle of the velocity vector to 
the orbit tangential during the thruster firing. Figure 5-18 shows the controls for this run. 
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Thrust vs Time 



Figure 5-18 

Excepting the delay with which the thrust starts its initial bum, the no theta control 
profiles are nearly identical to the original run. The initial thrust occurs with two distinct 
maximums and then approaches zero for the duration of the decay period. The final thrust 
is the characteristic fimal bum to match initial conditions. The angle alpha is nearly zero 
fi)r the duration of the thruster bums and wanders fi’eely between positive and negative pi 
when the thruster is off The cost fimction for the no theta run is 0.7945, which is only a 
1.4% reduction in optimality fi'om the original run. 

The final sinq)lification of the computer code involves the angle alpha. Each of the 
above cases shows that alpha hovers around zero whenever a thruster is fired. The logical 
procedure would then be to assume that the alpha is zero always. Physically this would 
mean that the thruster is constrained to fire in the direction of the velocity vector only. 
This assumption is commonly used in real world orbit reboost planning and execution. 

The removal of alpha allows the equations of motion to be dependent upon only one 
control. Every time alpha appeared in the equations of motion it was included within a 
trigonometric e3q>ression. Since alpha plays a bigger role while the thrusters are firing and 
the values were around zero, the trigonometric values of these regions were nearly one or 
zero depending on the trigonometric fimction. The programs notopt and notcrit were 
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exchanged for noaopt and noacrit, which both included substituting in all instances of 
alpha with the value zero. In all of the problems alpha was represented by the portion of 
the “aop” defined as “aop(5*n+i)”. Figure 5-19 shows the two lines of code in notcrit, 
which contain the alpha term. 

Ast = (aop(4*n+i)*cos(aop(5*n+i))-D(i))/(aop(3*n+i)*B); 

Ant = aop(4*n+i)*sin(aop(5*n+i))/(aop(3*n+i)*B); 

The expressions are replaced by 

Ast = (aop(4*n+i)*l)-D(i))/(aop(3*n+i)*B); 

Ant = 0; 

Figure 5-19 


The size of the “aop” vector is also reduced by “n” terms in the noaopt file. The first 
noticeable difference between the no alpha nm and the original run is the length of 
required computing time. The trimmed code and reduced matrices complete the 
optimization task in less than one half hour while the original code required over 2 hours 
to finish The results are shown below. Figure 5-20 shows the states of the system. 

Radius vs Time Velocity vs Time 














These graphs closely resemble the plots earlier shown for the original run. The cost 
functions are nearly identical with the no alpha case 0.03% more efficient that the original 
problem. The controls are shown in Figure 5-21. 


Thrust vs Time 



Alpha vs Time 
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Figure 5-21 

While the plot of alpha should be obvious, the plot of the thrust versus the time is almost 
identical to the thrust plot in Figure 5-4. 

The series of modifications to the program produced a code that minim ized the effects 
to the states yet modeled the system with a high degree of accuracy. Figure 5-22 shows a 
siunmary of the cost functions for each of the different modification steps. 


Method 

Cost 

Function 

Full Equations of Motion 

0.7837 

Fked Period 

0.8559 

No Theta 

0.7945 

No Alpha and No Theta 

0.7834 


Figure 5-22 
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The initial and final orbits of the optimal code are nearly circular with very small 
values for the flight path angle. The size of y(0) and y(T) are equal and are on the order of 
10"® radians. The values of the initial and final velocity also differ slightly from the 
reference value of one, usually in the sixth decimal place. The significance of these 
differences lies in the fact that the optimal initial conditions require the optimal initial orbit 
to be slightly elliptical. The no alpha program aUows the optimization code have a high 
degree of reproducibility with fewer errors or warnings given by MATLAB. Another note 
concerning the no alpha code is warranted. The plots displaying the controls of the 
program do show a curve for alpha, which is always zero for the no alpha code. The 
control history of the no alpha optimization code is used in the same first order propagator 
that the full optimization code uses. Displaying the control history of alpha reminds the 
viewer that the no alpha optimization code was used. The code also allows the tolerances 
for the constraints and states to be very small which in turn helps prevent the code fi'om 
settling in a local minimum without a full series of iterations. 

D. BAND FDCED SIMULATIONS 

The next type of trajectory to examine is the case where the spacecraft is bounded by 
an upper altitude limit. Following the rule of thumb previously established, a spacecraft 
subjected to this one constraint would follow a trajectory that would maintain the orbital 
altitude at the upper limit. This would be a thrust-drag cancellation profile and labeled 
earlier as the high FKT profile. Any flight below the upper limit would require more 
propellant. While appearing trivial the impact of the preceding statement is significant in 
that a satellite that possesses a maximum limit can not orbit the Earth with a trajectory that 
uses less propellant than a high FKT. To illustrate these statements the code was modified 
to make the initial altitude the upper limit as well. This modification forces the code to 
optimize be either following a trajectory that goes below the initial altitude or follow an 
FKT. Figure 5-23 shows the states for the results fi^om the FKT case. 
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Radius vs Time Velocity vs Time 





Figure 5-23 

While there may appear from the plot alone, that there are large oscillations in the 
radius during the beginning and end of the period, in actuality the peak of these 
oscillations has a magnitude of 1.0000001 or 10'’ deviation from the state. The 10'^ is the 
tolerance the optimal code uses in verifying the validity of its states. Thus neglecting these 
small numerical oscillations the radius follows an FKT profile. Figure 5-24 shows the 
controls for the FKT case. 

Once again neglecting the numerical inaccuracies at the very ends of the simulation the 


Thrust vs Time 



Figure 5-24 
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value of the thrust is a constant one. Since thrust is defined in terms of drag, a thrust 
value of one means that thrust equals drag. The satellite flies an FKT profile. 

However very few satellites today have the capability to follow a thrust-drag 
cancellation profile. Continuous or long duration bum thrusters are only beginning to be 
seen with the advent of electric propulsion systems. So most mission planners use the 
orbital band for planning the orbit maintenance schedule for a vehicle. To look at this type 
of mission the program first fixes the initial altitude and since the radius is periodic the 
final radius is also defined. The next parameter to examine is the period t. If x were fi'ee, 
the satellite would boost fi’om its lower altitude to its upper limit and then begin a high 
FKT for an infinitely long time. Since this defeats the purpose of the mission planner a tau 
needs to be chosen. The optimal code is given a fixed tau as in the earlier examples and 
fixes the initial altitude at an initial height. The only additional parameter is a value for the 
width of the orbital band, labeled “band” in the program noaopt. Figure 5-25 shows the 

g(i+6*n+5) = aop(i)-(l+band); 

Figure 5-25 


Radius vs Tame Velocity vs Tane 



Figure 5-26 

additional inequality constraint added to the noacrit program to boimd the upper altitude. 
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The first run was with the same period as before, 112.6 time units. The band was 
entered in as 0.001 DU or about 6.7 km. The initial altitude remained at 300 km and the 
generic spacecraft was the vehicle modeled. The no alpha code was used due to its 
numerical reliability and shorter computational time. The states for the first run are shown 
in Figure 5-26. 

The graph for the radius in the upper left shows the spacecraft follow a direct path to 
the high altitude limit, then begin an FKT trajectory, followed by decay to its original 
altitude. This combination of regimes can also be seen in the mass plot. Initially the rate 
of propellant use is large during the initial bum. The slope of the mass curve changes 
when then spacecraft reaches its upper altitude limit to a lower rate, consistent with a 
thmst-drag cancellation bum. The next change in the mass’ slope occurs when the slope 
becomes zero, during the decay phase when propellant is not used. The final change in the 
graph’s slope is due to the final bum required to match the initial conditions. Figure 5-27 
shows the controls for this case. 



Figure 5-27 

The thmst begins with a large initial bum with a gradual taper to a thmst 
approTdmately with the value of one. The thmst ends with the characteristic final bum. 

Numerous runs were made to examine to verify that the type of profile seen in the 
above example was consistent with varying orbital bandwidths and orbit maintenance 
periods. One should note that the constraint placed on the maximum radius was not 
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always used. If the band was sufficiently high or the period fairly short the satellite never 
made it to the point where it began an FKT profile. Instead the plots were similar to the 
plots shown in the fixed tau cases when the radius was free and unbounded. Whether or 
not a spacecraft’s radius was constrained depended on the period of the case. If the boost 
time and decay time filled the entire period without boosting the satellite to the constraint 
height an unbounded trajectory was the result. For completeness another example is 
shown. In this case the band is 0.005 distance units, just over 33 km. The period is 400 
time units. Figure 5-28 shows the states of this run. 


Radius vs Time Velocity vs Time 





Mass vs Trie 



Figure 5-28 


Figure 5-29 shows the controls of this example. Both figures resemble the state and 
control figures from the previous example. It is a reasonable conclusion that this is the 
shape of the optimal trajectory given the noted system constraints. 
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VI. ANALYSIS OF OPTIMALITY 


A. OPTIMALITY TEST CODES 


Previous sections have examined the optimal periodic problem, the numerical method 
and the implementation of the non-linear problem. This chapter is devoted to using the 
developed tools and examines their performance for varying different parameters. A 
measure of the performance is to analyze the optimality of the different solutions by 
comparing them to conventional maneuvers. A series of programs use the orbprop 
propagator to develop trajectories that simvdate maneuvers that are compared with the 
results of the optimal code. The first of these is called hoh and models the Hohmann 
transfer. This program is given an orbital bandwidth. It then calculates the change in 
velocity required to travel fi'om a circular orbit at the lower altitude to another circular 
oibit at the higher altitude. The Hohmann maneuver uses impulsive bums, meaning that 
the bums are instantaneous. To model this a velocity change is added directly to the 
vehicle when it reaches the bottom of the orbital band. At the same time the velocity is 
added, the mass used in the impulsive bum is subtracted fi’om the vehicle weight. Figure 
6-1 shows the lines of code which control the additions of velocity and subtractions of 
mass. 


First Impulsive Bum 

vfa = (l/xx(k,l))^.5; 

%Lower Circular Velocity 

vtxa = (2/xx(k,l)-l/atx)^.5; 

%Transfer Velocity at A 

dv = abs(vtxa-vfa); vn = xx(k,2)+dv; 


mf = xx(k,4)*exp(-dv/ve); 

%New Mass After Bum 

Second Impulsive Bum 

vtxb = (((2/(r+band))-(l/atx)))^.5; 

%Transfer Velocity at B 

vfb = (l/(r+band)yH).5; 

%Ifigher Circular Velocity 

dvb = abs(vfb-vtxb); 

mf = xp(z,4)*exp(-dvb/ve); 

%Fmal Mass after Transfer 

vn = dvb + xp(z,2); 

%Velocity addition 


Figure 6-1 
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The Hohmann program propagates the orbit over the given period. When the radius falls 
below the lower altitude limit, the Hohmann routine transfers the vehicle to the upper 
altitude. The mass and time of the Hohmann bums are recorded and then written to a data 
file for further use. 

The programs bchk and check are very similar and simulate a continuous bum type 
trajectoiy. The program check is given a maximum value of thmst. At the start of the 
simulation the vehicle is placed into a maximum tangential bum which continues until the 
vehicle breaks the top of the altitude band. The five physical states are then interpolated 
to give a specific time the spacecraft crossed the altitude upper limit. The next part of the 
program calculates the decay time from the top of the band to the lower limit. The times 
are added together to form a final period. The values for the mass and times for the one 
tangential bum method are recorded for latter comparison. The program bchk operates in 
much the same maimer with a minor exception. Here the program is given the entire 
period of interest beforehand. For example to compare the bchk trajectory versus the 
optimal code the same optimal time period x would be used. The program calculates the 
thmsting boost time and the decay time. The bum time and decay time are added and 



-20 0 20 40 60 80 100 120 140 

Time 

—♦—Optimal Code —■—Hohmann Transfer Check Program Band Check 

Figure 6-2 

subtracted from the given period. The time remaining is fiJled with a high FKT at the 
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upper altitude limit. The trajectory resembles the band-limited trajectories included in the 
last chapter. Here also the values of the masses and times are recorded for comparison. 

Figure 6-2 shows the differences in the mass plots of the four different trajectory 
programs. The two check programs require considerably more fuel than the Hohmann 
and Optimal trajectories. The codes are meant to mimic a real world continuous thrust to 
the upper altitude limit. The trajectories are very costly in terms of propellant due to the 
fact of the bum required to circularize the orbit to the upper altitude limit. There are two 
reasons that this needs to be accomplished. The first is that the upper circular orbit is the 
same final orbit after the Hohmann transfer. The second reason is that the angle gamma is 
a positive number, meaning that the velocity vector is above the local horizontal. After 
the thmster is turned off the spacecraft would travel beyond the upper altitude constraint 
in an elliptical oibh. Figure 6-3 is a diagram showing the necessary change of velocity 
required to circularize the final orbit at the upper altitude limit. 



For low values of thmst, where T = 1 to 5 normalized ffirust units, the required Av is fairly 
small. When thrust increases, y increases lengthening the required velocity change vector. 
Thmster sizes that approximate instantaneous bums require a large second bum to 
drcularize the final orbit. The thmster size used in figure 6-2 for the check and bchk 
programs was set at 20 normalized thmst units. All four programs will be examined later 
to measure each trajectory’s optimality. 

B. VARIATION OF PARAMETERS 

The tme test for the code is its ability to vary the inilial guesses, control boimdaries, 
numerical method parameters and deliver reliable and accurate results. This is especially 
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crucial for optimal control problems that are very sensitive to initial values and 
effectiveness of the controls. In the previous chapters different values and equations were 
tested to achieve good results. Throughout these examples, thrust has been constrained 
to never go above a certain value. In this case thrust has been constrained with a value of 
5 normalized units. Physically the value of thrust changes with the location and physical 
attributes of the spacecraft. For example for the generic spacecraft at 300 km altitude, 
one normalized unit of thrust equals 0.7 N while the same normalized unit for the ISS at 
the same altitude is 2.9N. The engine plaimed for the ISS is 1 ION or just over 37 
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Figure 5-19 

normalized units at the reference altitude. Therefore it is beneficial to study the effects of 
increased thrust constraints and eventually the unconstrained thrust case. 

Recalling Figure 5-19, the controls of the generic satellite given a fixed tau of 112.6 
TU, the plot for the thrust shows the smooth continuous thrust history with the initial 
thrust defined by two characteristic peaks with a maximums around 5.0. Note that the 
reason the graphs appear to shghtly break the upper thrust constraint are due to the fact 
that the plots represent an interpolation of the values taken fi-om the values of the controls 
determined by the optimization code at the LGL points. In attempting to smooth the 
curve the interpolation routine will go above the stated constraints. However the shape of 
the curve is still a good approximation of the behavior of the system. Wfith this thrust 
constraint the cost function for the given period was 0.7834. The propellant consumed 
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during the simulation was approximately 0.0057 normalized units or in physical terms just 
under 17 kgs. 

Thrust was then increased to a value of 10 thrust units. This must be done in two 
places within the computer code. First in the noaopt code the initial thrust vector was set 
to 10 units. In the noacrit code the inequality constraint containing “aop(4*n-i-i)” was 
given a maximum limit of 10. Figure 6-4 shows the control history of the new simulation. 
This figure shows another smooth function for the thrust but with a single peaked initial 


Thrust vs Time 



Figure 6-4 

bum. The final bums at the end of the simulation are about the same for both cases and 
do not reach the maximum thrust value. The cost function for this case is 0.7688, a 2% 
reduction in fuel consumption. In terms of physical mass of fuel saved by having a higher 
limit of thmst, it represents 0.3 kg of savings. This may not seem like a substantial 
amount but recall that the period for this case is slightly longer than a single day. 
TnfTftasin g the motor size fi-om 3.5 N to 7 N would save the spacecraft approximately 100 
kg of propellant needed for orbit maintenance over the course of a year. The reason for 
the fuel savings is that the larger motor is able to boost the satellite to a slightly higher 
orbit with the period constraint. The 5 thmst unit case boosted the satellite to a 
normalized altitude of 1.0034 DU while the 10 thmst unit case boosted it to 1.0036 DU. 
This change in altitude difference was enough to reduce the cost function ever so slightly. 
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The shorter bum time to reach the higher altitude allows longer decay times. The basic 
shape of the trajectory is very similar to the earlier case. Figure 6-5 shows the states for 
the 10 thmst unit example. 


Radius vs Ttne Velocity vs Time 



Figure 6-5 

The states have the same shape as those of the lower thrust case, so the addition of thmst 
does not alter the trajectories significantly. 

Since the thmst history was constrained at the maximum thmst limit during the last 
simulation, another run is necessary with higher thmst levels. The next simulation used a 
maximum available thmst of 25 thmst units (14 N). The ultimate goal of raising the thmst 
constraint is to find a case where the maximum thmst does not reach the constraint; in this 
sense, the thmst would be xmbounded. An easy option would be to set the thmst limit at a 
very high number in the hundreds or thousand. Unfortunately the code sees the large 
range of available thmsts as an added complexity for the system, leading to non-smooth 
and non-optimal solutions. Thus an incremental approach to raising the altitude limit is 
desired. The controls for this run are included in figure 6-6. 

The controls fi’om the T=10 and T=25 examples are very similar. The most significant 
point of figure 6-6 is that maximum thmst is 14.3 thmst units. This is the first case where 
the thmst did not reach the maximum constraint with an optimal trajectory. The cost 
function of this mn is 0.7620. This cost function is a 3% reduction over the original thmst 
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constrained case. To confirm that this solution is not a numerical accident another case 
with a different thrust constraint (T=20) was run. It is important to determine that the 
optimal thrust profile is not a function of the location of the thrust constraint or initial 


Thrust vs Time 



Figure 6-6 

guesses for the thrust history. The cost function of this second run was also 0.7620. All 
of die states were the same with minor fluctuations in the thrust history in the 5^ decimal 
place. 

Figure 6-7 shows the controls for the T=20 nm. 


Thrust vs Time 



Figure 6-7 
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The thrust functions are identical between the two runs. The states are included in 


Figure 6-8. 


RadiiB vs Time NAslodty vs Time 




The states shown in Figure 6-8 represent the optimal solution given the generic spacecraft, 
starting at the reference altitude and given a certain period. This procedure can be 
repeated for any variations of the given information. However, the states shown in Figure 
6-8 represent the optimal paths for this unique situation. To develop other situation only 
requires changing the conditions under which noaopt and its associated codes run. For 
example a different set of curves can be determined if the satellite was to have a two and a 
half day orbit maintenance period. 

The next step is to check the results of the optimal code against the Hohmann 
trajectories given the same system constraints. One of the tenets fi’om the beginning of the 
study was to identify a trajectory that was as good or slightly better than an ideal 
Hohmann trajectory. Figure 6-9 shows the mass consumption of various profiles including 
the optimal code run with T^ax = 20, the Hohmann transfer and several FKTs. 
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The FKT paths represent, from the top to bottom in Figure 6-9, the low-band, mid-band. 


Propelart Comparison 



Figure 6-9 

and high-band trajectories. As noted earlier the high-band FKT will always be the 
trajectory with the lowest fuel consumption but the figure assumes that the satellite begins 
the simulation at the top of the band, neglecting the boost maneuver needed to reach the 
top of the band. The straight stair-step like solid line represents the mass consumption of 
the Hohmann transfer. The line with the plus signs represents the mass consumption of 
the optimal trajectory. Notice that optimal code shows lower propellant use. In order to 
highlight these differences the graphing periods are increased to show savings over longer 
periods of time. Figure 6-10 shows the extended mass consumption comparison for the 
constrained thrust example. 
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Mass Consumption Comparison 



Time 

-Optimal Code -Hohmann Transfer 

-Check Program -Band Check 

Figure 6-10 

The uppermost line shows the propellant consumption for the two check programs. In 
this instance both programs show the same trajectory. This is because the time to decay 
from the top of the band is greater than the given orbit maintenance period. When the 
decay time is subtracted from the given period the result is a negative number. The bchk 
code considers this value zero; in other words the time that the satellite follows a high- 
band FKT is zero. The next two lines are very close together but a difference is 
discemable. The higher of the two lower lines is the propellant use of the Hohmann 
trajectory. The optimal trajectory has the lowest rate of fuel use of the four tr^'ectories. 
In order to highlight the differences between the optimal and Hohmann trajectories a 
trendline is extended from each set of points. Figure 6-11 shows the trendlines for a 
period of approximately 150 days. 
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Mass Consumption Comparison 



Optimal Code -Hohmann Transfer 

Linear (Hohmann Transfer)-Linear (Optimal Code) 


Figure 6-11 

This graph shows that one of the goals of this project has been met. There exists a 
trajectory that is more fuel-efficient than a Hohmann transfer. While the difference 
between the two is small one must remember that the Hohmaim transfer was modeled as 
an ideal maneuver. The bums are instantaneous with infinite thmst. This assumption 
neglects losses, mainly drag loss, which the optimal control code considers automatically. 

The easiest way to show the effect of drag loss is to start with Newton’s equation of 
motion with the assumption that mass changes are small. 

F^=ma (76) 

Assuming constant acceleration equation (76) becomes; 


m 


= Av 


F„ = T-D 


(77) 

(7S) 


Thus the required Av can be related to the total impulse through equation (77). Recalling 
a form of the rocket equation. 
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mp = m^ 


-iiv 


1 ~ ^ 


(79) 


For a real motor: 


/w =• 


Assuming constant thrust: 


m = - 


So^sp 


At 


(80) 


(81) 


Here At is the bum time of the thruster. Combining equations (80) and (81) gives an 
expression for total impulse in terms of propellant mass. 


FAt = 




So^sp 


(82) 


Assuming that drag works in opposition to thrust, equation (82) can be rewritten. 


TAt-DAt = 




So^sp 


(83) 


The DAt term from equation (83) is the portion of the total impulse lost to drag. Since the 
total impulse must remain the same for a given change in velocity, the drag impulse must 
be made up by the thrust impulse. This either requires a larger thruster or a longer bum 
time. For example, assume that a total impulse of200 units is needed. With a thruster 
with a normalized thmst of 2, or twice the force of drag, the spacecraft would need a 
thruster firing of 100 time units to accomplish the maneuver if there were no drag. With 
drag, the net force is equal to one unit and requires a bum time of two himdred time units. 
The extra time required for the maneuver determines the amoimt of propellant required. 
The addition of drag uses twice as much propellant, or in other words the drag loss of this 
example is approximately 50%. With larger thrusters the drag loss is decreased. However 
even for a motor that has a rating of 15 normalized units, the optimal thruster example 
determined previously, nearly 7% more propellant is required by considering drag. Figure 
6-12 modifies the Hohmann propellant consumption of figiare 6-11 by adding a 6.67% 
increase in the amount of fuel required. 
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Mass Consumption Comparison 



-Optimal Code - 

-Hohmann Transfer 

-Linear (Hohmann Transfer)- 

-Linear (Optimal Code) 


Figure 6-12 

Figure 6-13 is included in order to show numerical values for the fuel savings of the 
optimal code. 



Mass 

Expended 

Percent 

Savings 

Launch 

Cost 

Savings 

optimal Code 

644.3 kg 

- 

- 

Ideal Hohmann 

667.0 kg 

3.5% 

-$250,000 

Hohmann with 

Drag Loss 

711.5 kg 

10.4% 

-$740,000 


Figure 6-13 


The sampling point for Figure 6-13 was at approximately 45 days of orbit maintenance 
operations. The thruster size was set at the maximum of the optimal code, 15 normalized 
thrust units. The 3.5% savings the optimal trajectory has over the impulsive Hohmann by 
itself is a significant improvement over the life of the satellite. When drag effects are taken 
into consideration the savings improve dramatically. The general rule remains that the 
smaller the orbital maintenance thruster the greater the savings of the optimal code 






compared to the finite-bum Hohmann transfer. The above example shows that 10.4% 
savings in propellant can be achieved when the size of the motor equals 15 thmst units 
(10.5 N). The efiBciency of the Hohmann transfer can be improved by increasing the size 
of the orbit transfer engine, which decreases the bum time. However even at the 
maximum limit, an infinite thmster, the Hohmann transfer is still 3.5% less eflBcient. 

When monetary terms are used to describe the potential savings the differences 
become more evident. Looking at Figure 6-13, the savings between the Hohmann and the 
optimal code is 23.3 kg of fuel. Using an average launch cost to LEO this comes to 
approximately $256,000. When drag losses are considered, the difference increases to 
67.2 kg, or in terms of dollars $740,000. The generic satellite used in the example is a 
high drag vehicle that operates at a fairly low altitude. However, sa\nngs occur with 
other spacecraft as weU. For example, the mass/area ratio was changed to model the ISS 
during a late stage in constmction. Figure 6-14 shows the physical characteristics of the 
ISS example. 


ISS Characteristics 

Orbital Altitude 

300 km 

Mass 

408,420 kg 

Area 

2200 m^ 

Drag Coefficient 

2.35 


Figure 6-14 


One rather large difference between the future ISS operations and the optimal code is 
that the length of time between rendezvous for the ISS is much longer than the period 
entered into the optimal code. The main reason for this is that long periods create some 
numerical instability for the optimal code without extensive gradient information. The 
example though does highlight some of the benefits of an optimal trajectory for shorter 
time periods versus the mission planned with Hohmann type transfers. All that was 
changed fi^om the generic satellite to the space station example were the constants for 
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mass and area, included in the beginnmg of the noaopt program. A maximum thrust of 20 
was placed as the initial guess for the code, which for the ISS is approximately 57.7 N. 
The period of the optimal code was fixed at 112.6, the same period as previous examples. 


Radius vs Time 




Velocity vs Time 



Figure 6-15 

Figure 6-15 shows the states for this example. 

The plots for the states are nearly identical in shape to the states of the generic satellite. 
The difference lies in the values of the states. The ISS is much more masrive than the 
generic vehicle and has a normalized ballistic coefiBcient two order of magnitudes larger. 
Figure 6-16 shows the controls for this run. 


Thrust vs Time 



Figure 6-16 
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Again Figure 6-16 does not possess any large deviations from earlier plots. One 
interesting point to notice is that the maximum thrust for this case is just over 10 units (or 
about 30 N). For this case and period the optimal thruster is much smaller than the 
thruster currently being planned for the ISS, by approximately a factor of four. Figure 6- 
17 compares the Hohmann transfer fuel consumption versus the optimal trajectory. 

xIO"* Propellant Comparison 



Figure 6-17 


The optimal code is the hash-marked line that increases smoothly as the vehicle undergoes 
its transfer bum and then is constant during the following decay. The figure also shows 
two different Hohmann transfer plots. The solid stair-step line is the mass consumption of 
the Hohmann transfer when the transfer uses the same size orbital band as the optimal 
code. The dashed stair step line is a Hohmann transfer where the period of the transfer is 
the same as the period of the optimal code. The straight lines are the three different FKT 
profiles. Since the band in this case is very narrow (about 1.5 km) there is little difference 
between the high-band to low-band FKT. The fixed tau Hohmann transfer nearly matches 
the optimal fuel efficiency until the final circularizing bum at the end of the period. At 
first glance Figure 6-17 suggests that the fixed band Hohmann transfer uses less propellant 
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than the optimal code. However, the period of the fixed band Hohmann transfer is shorter 
than the optimal trajectory. In order to examine the difference between the two orbit 
transfers the values of each are linearly extrapolated over a period of 150 days. This 
highlights long-term divergence between the two methods. 



Fi^re 6-18 


Figure 6-18 shows that again the optimal transfer is a more efficient maneuver than an 
impulsive Hohmann transfer. For the maximum thruster, of 10 units, drag losses would 
increase the consumption of propellant by the Hohmann transfer by over 10%. In terms of 
actual propellant, a sample at the 50-day mark would show that the optimal trajectory 
used 4252.9 kg of fiiel for orbit maintenance while the Hohmann transfer would use 
4337.2 kg. The savings in this case are approximately 2% over the ideal Hohmann, drag 
loss may increase that number to 13%. The optimal code would save about 600 kg per 
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year over the ideal Hohmaim transfer, a savings of 6.6 million dollar per year in laimch 
weight costs alone. 

This example shows the usefulness of the code as well as the ease with which different 
spacecraft can be modeled. While a lower drag vehicle like the ISS has slightly less 
savings over the ideal Hohmann than a higher drag vehicle, the savings are present and 
significant over a long period of time. Further analysis is necessary to vary the parameters 
and understand the “physics” of the optimal trajectory. 


80 


VII. CONCLUSION 


The work succeeded in accomplishing the stated goals. Quantitatively, the rule of 
thumb is that the optimal trajectory is slightly better than the impulsive Hohmann (~ 3%) 
and significantly better than the finite-bum Hohmann (up to ~ 10%). The savings are at 
least equal to drag loss “elimination.” Beyond the physical numbers, the cracial part of the 
thesis was to explore the orbit maintenance problem in depth and use optimal periodic 
control (OPC) theory as a tool to locate the fuel-optimal solution. The many simulations 
that were described were a very small portion of the learning and experimentation that 
went into attacking this process. The work began with general optimal control theory, 
later specialized by periodicity, the equations of motion for two-dimensional orbital travel, 
and a good numerical method. Examining the results of the preliminary simulations reveal 
some numerical instabilities in the optimization code which led to further examination of 
the equations in order to improve the reproducibility and accuracy of the results. With this 
accomplished, the optimal results could then be compared with various other trajectories 
in order to measure the optimal trajectory’s fuel efficiency. The different controls and 
boundary conditions were modified in order to search for the most fuel-optimal trajectory. 

The amoimt of iirformation fi'om these simulations is enormous. The optimal 
tr^ectory would rqrresent a large departure fi'om current orbit maintenance procedures. 
The first difference would manifest itself in the design process. Mission planners 
determine the periodicity of the orbit maintenance routine by either fixing an orbital band 
in which the satellite operates or by specifying specific times the spacecraft is to be at the 
lower altitude limit. With this information, the operating altitude, and the satellite’s 
physical characteristics the optimal code can be run with varying thrust levels. By 
examining the control histories a maximum thrust can be determined. This is an absolute 
limit for the given conditions. Access performance in the engine would not increase fuel 
efficiency. This limit gives designers an upper limit on the size of the needed thruster, 
creating possible savings of platform cost and weight. 

Then there is the trajectory itself A true departure fi'om the Hohmann transfer, the 
optimal code calls for a smooth continuous throttle bum that transfers the vehide to the 
top of the orbital band. The only requirement for the bum would be that since the burp 
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occurs over a larger period of time, it should be performed as an autonomous maneuver or 
be continually linked to the control of a ground station. Most current engines are 
normally bang-bang only. This trajectory ^ves enough justification by its savings to 
develop throttlable motors or develop chattering techniques which can model bang-bang 
engines as a varying thrust engine. 

But with all good research there are new questions that can be addressed and studied. 
First the complexity of the atmospheric model can be increased. Diurnal effects may be 
able to increase the savings of trajectories. One of the strengths of this method is that 
modifying the atmospheric model only requires modification of the density function that is 
independent fi’om the optimization code. Another area of future research is identifying 
more robust optimization codes and routines. A good portion of this work was spent in 
determining ways in which to increase the numerical accuracy of the results while also 
reducing the computer run time. New methods for numerical approximations may be 
highly productive and more accurate interpolation schemes may have a large benefit. The 
code is not limited to Barth bound operations. Simply by chan^g the constants 
throughout the programs and the density model, this method could be applied to any other 
planetary body. 

The biggest strength of the program is its flexibility and applicability to (fififerqnt 
models. The normalization method made it very easy to modify the spacecraft 
characteristics or the size of the thrusters. In terms of physical characteristics, spacecraft 
are modeled by their mass to cross-sectional area ratios. This is what makes the ISS the 
same as the Hubble Space Telescope but different firom an inflatable. The id)ility to run 
the code for different situations is an important aspect of the project’s success. 

Optimal control has long broken preconceptions about certain topics. The problem of 
orbit mainteiumce may be the latest example. In every primary orbital mechanics textbook 
there is a section on the ^lohmaim transfer that states that it is the most efficient transfer 
between circular dibits, which it is for many cases. But when complexities are added to 
the system, such as the addition of atmospheric drag, these statements do not hold true. It 
is the goal of every research project to question all the answers and push the boundaries of 
learning and understanding out a little further. 
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Appendix A - Density.m 


APPENDICES 


% Program density.m 
% Karl Jensen - May 13,1998 

function rho = density(x); 
rt = x* (6378+300); 

if rt >= 6753 

rref = 6778; h=58.2; 
rhoref = 2.62e-12; 
elseif rt >= 6703 

rref= 6728; h =54.8; 
rhoref = 6.66e-12; 
elseif rt >= 6653 

rref = 6678; h=50.3; 
rhoref = 1.87e-ll; 
elseif rt >= 6603 

' rref = 6628; h =44.8; 
rhoref=5.97e-ll; 
else rref = 6578; h= 37.5; 
rhoref = 2.41e-10; 

end; 

rhon = rhorePexp(-l *(rt-rreQ/h); 
rho = rhon/1.87e-ll; 

Appendix B - Orbprop.ni (First Order Propagator) 

% First Order Propagator 
% Karl Jensen - May 1,1998 

function ddot = orbprop(t,x) 

global T alpha tau n band xt 
global xt 

% First Order Differential Equation Solver 
% States 

% x(l) = r 

% x(2) = v 

% x(3) = gamma 

% x(4) = mass 

% x(5) = theta 

m0 = 408420; a=2200; 

%m0 = 3000;a=1000; 



m0=3000;a=500; 

isp = 300; gn = 9.81; cd = 2.35; 

rok = 6378+300; rhoO = 1.87*10^-11); rho =1; 

vO = sqrt((3.986005*10^5)/rok)*1000; rO = rok * 1000; 

t 

B = 2 * mO / (rO*rhoO*cd*a); 
ve = isp*gn/vO; 


rho = density(x(l)); 

g=l/(x(ir2); 

D = rho*(x(2)^2); 
if (n>5) 

% I assume lower limit of the run = 0; 
taut = (2*t-tau)/tau; 
t 

Tp = interpl(xt,T,taut,'spline'); 
ifTp<0 

Tp = 0; 

end; 



alphap = interpl(xt,alpha,taut,'spline'); 

end; 

if (n=l) 

Tp = D; alphap = 0; 

end; 

if(n=2) 

Tp = 0; alphap = 0; 

end; 

if(n=3) 

Tp = 20; alphap = 0; 

end; 



As = (Tp*cos(alphap)-D)/(x(4)*B); 

An = Tp*sin(aiphapy(x(4)*B); 

ddot(l) = x(2)*sin(x(3)); 

ddot(2) = -g*sin(x(3))+As; % 5 Equations of Motion 

ddot(3) = (x(2)^2/x(l)-g)*(cos(x(3))/x(2))+An/x(2); 
ddot(4) = -abs(Tp)/(ve*B); 
ddot(5) = x(2)*(cos(x(3))/x(l)); 
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APPENDIX C - HOH.M (HOHMANN TRANSFERS) 

% Drag Code - Karl Jensen - Orbital Mechanics November 26,1996 
global T 

global band n tau ve 
mu=l; 

%band = 1.000147536-1; 
m0 = 408420; a=2200; 

%m0 = 3000; a =500; 

%m0=3000;a=10000; 

isp = 300; gn = 9.81; cd = 2.35; 

rok = 6378+300; rhoO = 1.87*10^-11); rho =l;radref =rok; 

. vO = sqrt((3.986005*10^5)/rok)*1000; rO = rok * 1000; 


ve = isp*gn/v0; 


% t,x are for Constant Thrust - Lower FKT 
% tt, XX are for Free Fall 
% tl, xl are mid band FKT values 

1 = 1 ; 

v=l; 

gm=0; 
m=l; 
theta = 0; 

p = 2*pi*(r^3/mu)^.5; 
tp = tau/(2*pi)*p; 

% Constant Thrust 
xo=[r V gm m theta]; 
to = 0; 
tf=tp; 
tol= l.e-9; 
n= 1 

[t,x] = ode45('orbprop',to,tf,xo,tol,0); 

figure(l); 
orient tall; 

subplot(2,2,1 ),plot(t/p,x(:, 1 ),'b')pdabel('Orbits');title('Radius versus Orbit');ylabel('Radius 

(m)'); 

subplot(2,2,2),plot(t/p,x(;,2),T5');xlabel('C)rbits');title('Velocity versus 
Oibit');ylabel('Velocity (m/s)'); 
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subplot(2,2,3),plot(t/p,x(^3)*360/(2*pi),'b');xlabel('Orbits');title('Gamma versusOrbit'); 
ylabel('gajnma (deg)'); 

subplot(2,2,4),plot(t/p,x(:,4),'b')pdabel('Orbits');title('Mass versus Orbit');ylabel('Mass 
(kg)'); 

% Free Fall 
n = 2 
sp = 5*p; 

[tt,xx] = ode45('orbprop',to,sp,xo,tol,0); 

temp = xx;temp2 = tt; 

figure(2); 
orient tall; 

subplot(2,2,l),plot(tt/p,xx(:,l),'b');xlabel('C)rbits');title('Radius versus Orbit');ylabel('Radius 

(m)'); 

subplot(2,2,2),plot(tt/p,xx(:,2),'b');xlabel('Orbits');title('Velocity versus Orbit'); 
ylabel('Velocity (m/s)'); 

subplot(2,2,3),plot(tt/p,xx(:,3)*360/(2*pi),'b');xlabel('Orbits’);title(’Gamma versus Orbit'); 
ylabel('gamma (deg)'); 

subplot(2,2,4),plot(tt/p,xx(;,4),'b')pdabel('Orbits');title('Mass versus Orbit');ylabel('Mass 
(kg)'); 

% Middle Band FKT 
n=l; 

mid = r + band/2; Iv = (mu/niid)^.5; 
xo = [mid Iv gm m theta] 

[tl,xl] = ode45('orbprop',to,tf,xo,tol,0); 

% High Band FKT 
n=l; 

higher = r + band; hv = (mu/higher)^.5; 
xo = [higher hv gm m theta] 

[th,xh] = ode45('orbprop',to,tf,xo,tol,0); 


%Part Three - Hohmann Runs 


flag= l;n = 2;c=l; 


while flag 

k = c; 

while (xx(k,l)>(r)) & flag 
k = k+l; 
if k >= size(tt,l) 
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flag = 0; 

end; 

end; 

if (c 1) & flag 

buravec = xx(k-l:k,;); burnt = tt(k-l;k); 
int = mterpl(bumvec(;,l),bumt,l); 
invec = interpl(bumvec(:,l),bumvec,l); 
xx(k,:) = delvee(invec); tt(k) = int; 


end; 

if flag 

dr = pi*(cd*a/(m0*xx(k,4)))*rho0*(radref''1000*xx(k,l))'^2; 

dm = dr / (radrePlGOO); 

atx = (l+r+band+0)/2; 

vfa = (l/xx(k,l))^.5; 

vtxa= (2/xx(k;l)-l/atx)^0.5; 

dv = abs(vtxa-vfa); vn = xx(k,2)+dv; 

mf = xx(k,4)*exp(-dv/ve); 

n=2; txoib = [xx(k,l) vn xx(k,3) mf theta]; 

tf=tt(k) + pi*(atx^3)^.5; 

[tx,xp] = ode45('orbprop',tt(k),tf,txorb,tol,0); 
z = size(tx,l); 

vtxb = (((2/(r+band))-(l/atx)))^.5; 

vfb = (l/(r+band))^.5; 

dvb = abs(vfb-vtxb); 

mf = xp(z,4)*exp(-dvb/ve); 

vn = dvb + xp(z,2); 

to = tf; tf = tp; 

orb = [xp(z,l) vn xp(z,3) mf xp(z,5)]; 
ovr = 6*(band+0.005)/(dm); 
if (tf-to) > ovr 

tf = to + ovr; 

end; 

[tq,xq] = ode45('orbprop',to,tf,oib,tol,0); 
fif = [5 1]; 
if flr= size(xq) 
xq=xq'; 

end; 

St = [tt(l;k); tx; tq]; 
sx=[xx(l;k,;); xp;xq]; 
tt = st; 
xx = sx; 
c = k+z; 

end; 

end 
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figure(3); 

xx=temp;tt=temp2; 
orient tall; 

subplot(2,2, l),plot(st/p,sx(;, 1 ),V); 

xlabel('Orbits');title('Radius versus Orbit');yIabel('Ra(iius (m)'); 
subplot(2,2,2Xplot(st/p,sx(:,2),'b');xlabel('Orbits'); 
title('Velocity versus Orbit');ylabel('Velocity (m/s)'); 
subplot(2,2,3),plot(st/p,sx(;,3)*360/(2*pi),'b'); 
xlabel('Orbits’);title('Gamma versusOrbit');ylabel('gamma (deg)'); 
subplot(2,2,4),plot(st/p,sx(;,4),'b')pdabel('Oibits'); 
titleCMass versus Orbit');ylabel('Mass (kg)'); 
sxl =sx;stl=st; 


figure(4); 
orient tall; 

subplot(2,2,1 ),plot(st/p,sx(:, 1 ),'b'); 

xlabel('Orbits');title('Radius versus Oibit');ylabel('Radius (m)'); 
subplot(2,2,2),plot(st/p,sx(;,2),'b');xlabel('Oibits'); 
title('Velocity versus Orbit');ylabel('Velocity (m/s)'); 
subplot(2,2,3),plot(st/p,sx(;,3)*360/(2*pi),'b'); 
xlabel('Orbits');title('Ganima versusOrbit');ylabel('gamma (deg)'); 
subplot(2,2,4Xplot(st/p,sx(:,4),'b')pdabel('Orbits'); 
titleCMass versus Oibit');ylabel(Mass (kg)'); 


figure(5) 

plot(st 1 /p,m-sx 1 (: ,4),'b',t/p,m-x(: ,4),'m—',tl/p,m-xl(; ,4),'m: ',th/p,m-xh(: ,4),'m- 
.');title('PropeIlant Comparison'); 
xlabel('Orbits');ylabeI('Normalized Mass'); 
hold; 

%plot(stl/p,m-sxl(;,4),'r'); 
hold oflF; 

k=size(xl(:,4)); mid = m-xl(k,4) 
k=size(xh(:,4)); high = m-xh(k,4) 
k=size(x(:,4)); low = m-x(k,4) 

hohfuel = m-sxl(:,4);lowfuel = m-x(:,4); 
midfuel = m-xl(:,4);highfuel = m-xh(:,4); 
hohtime = stl; 
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Appendix D - Orboptm 


% This routine calculates the minimum of the cost function given by 

% the function crit. 

clear; 

tic; 

global n 

global Dn X w 

global B ve rho aop aref 

n=24; % The number of lobatto points. 

%m0 = 408420; a=2200; 
mO = 3000; a =500; 

%m0=3000;a=10000; 

isp = 300; gn = 9.81; cd = 2.35; 

rok = 6378+300; rhoO = 1.87*1(KX-11); rho =1; 

vO = sqrt((3.986005*10'^5)/rok)*1000; rO = rok * 1000; 

B = 2 * mO / (r0*rho0*a*cd); 


ve = isp*gn/v0; 


Tguess = 5.0; 

[Dn,x,w]=di^(n); 

aop(l;n)= ones(n,l); 
aop(n+l:2*n) = ones(n,l); 
aop(2*n+l:3*n)= l.l*ones(n,l); 
aop(3*n+l:4*n) = ones(n,l); 
aop(4*n+l:5*n) = Tguess*ones(n,l); 
aop(5*n+l;6*n) = 0*ones(n,l); 
aop(6*n+l:7*n) = zeros(n,l); 
aop(7*n+l)=112; 
taup=aop(7*n+l); 
options(l 3)=5 *n+6; 
options(14) = 400*(7*n+l); 
options(4)=l 0^-4);%5 
options(3)=l 0^(-5);%5 
options(2)=l 0^-5);%5 
options(l)=l; 

%options(18)=l; 
options(16) = 10^-6);%7 
vub=[]; 
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vlb=[ ]; 

aop=constr('dorbcrit',aop,options,vlb,vub); 

% a is for the state and b is for the control. 

rp=aop(l:n); 

vp = aop(n+l:2*n); 

gammap = aop(2*n+l:3*n); 

massp = aop(3*n+l;4*n); 

Tp = aop(4*n+l:5*n); 
alphap = aop(5*n+l:6*n); 
thetap = aop(6*n+l:7*n); 
tau = aop(7*n+l); 

toe; 

tic; 

orbres; 

oihconv; 

toe; 

Appendix E - Orberitm 

function [costfii,g]=dorbcrit(aop); 

%% This function calculates the cost function that is to be minimized 
%% and the state constraints. 

global nx wDn; 
global B ve rho aref; 

% global costfii; 

% Set up the cost function, 
for i=l:n 
fii(i)=aop(4*n+i); 
end; 

costfii= l/(2)*sum(w.*fii'); 
tau = aop(7*n+l); 


% Set up the state constraints, 
for i=l:n 

rt = aop(i)*(6378+300); 
if rt >= 6728 
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rTef= 6778;h=58.2; 
rhoref = 2.62e-12; 
elseif rt >= 6628 

iTef= 6678; h =50.3; 
rhoref = 1.87e-ll; 
else rTef= 6528; h= 37.5; 
rhoref = 2.41e-10; 

end; 

rhon = rhoref'exp(-l*(rt-rreQ/h); 
rho = rhon/1.87e-ll; 

G(i) = l/(aop(i)^2); 

D(i) = rho*(aop(n+i)^2); 

Ast = (aop(4*n+i)*cos(aop(5*n+i))-D(i))/(aop(3*n+i)*B); 

Ant = aop(4*n-H)*sin(aop(5*n-H))/(aop(3*n-i-i)*B); 

V(i) = ve; 

% Theta 

g(i+4*n) = (2/tau)*sum(Dn(i,:).*aop(6*n+l ;7*n))-(aop(n+i)/aop(i)*cos(aop(2*n-Hi))); 
%Mass 

g(i+3 *n)= (2/tau)*sum(Dn(i,:).*aop(3 *n+l :4*n))+(aop(4*n+i)/(ve*B)); 

% Gamma 

g(i+2*n)= (2/tau)*sum(Dn(i,:).*aop(2*n+l :3 *n)); 

g(i+2*n^g(i+2*n)-(((aop(n+i)^2)/aop(i))-G(i))*(cos(aop(2*n+i))/aop(i+n)); 
g(i+2*n^ g(i+2*n)- .^it/aop(n4i); 

% Radius 

g(i)= (2/tau)*sum(Dn(i,;).*aop(l:n))-aop(n+i)*sin(aop(2*n+i)); 

% Velocity 

g(i+n)= (2/tau)*sum(Dn(i,:).*aop(n+l:2*n))+(G(i)*sin(aop(2*n+i)))-Ast; 


g(i+5*n+6) = - aop(4*n+i); 
g(8*n+6+i) = aop(4*n+i)-5.0; 

% g(i+9*n+6) = l-aop(i); 
%g(i+7*n+5) = aop(i)-l; 
g(i+6*n+6) = aop(5*n+i)-pi; 
g(i+7*n+6) = -pi-aop(5*n+i); 
end; 

g(5*n+l) = aop(l)-aop(n); 
g(5*n+2) = aop(2*n+l)-aop(3*n); 
g(5*n+3) = aop(6*n+l); 
g(5*n+4) = aop(n+l)-aop(2*n); 
g(5*n+5) = aop(3*n+l)-l; 
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g(5*n+6) = aop(l)-l; 
g=g'; 

Appendix F - Orbconv.m 


% Conversion from AOPs to Functions of Time 

global aop n 
global T alpha tau 
global band xt 


rO = aop(l); vO = aop(n+l); gmO = aop(2*n+l); 
mO = aop(3*n+l); thO = aop(6*n+l); 

T = aop(4*n+l:5*n); 
alpha = aop(5*n+l :6*n); 

tau = aop(7*n+l); 
xt=x; 

xo=[rO vO gmO mO thO]; 
to = 0; 
tf = tau; 
tol= l.e-9; 

[t,x] = ode45('orbprop’,to,tf,xo,tol,l); 
nl = size(t); 

rr(l) = aop(l); rr(nl(l,l)) = aop(n); 
vr(l) = aop(l+n); vr(nl(l,l)) = aop(2*n); 
gmr(l) = aop(l+n*2); gmr(nl(l,l)) = aop(3*n); 
massr(l) = aop(l+3*n); massr(nl(l,l)) = aop(4*n); 
thrustr(l) = aop(l+n*4); thrustr(nl(l,l)) = aop(5*n); 
alphar(l) = aop(l+5*n); alphar(n 1(1,1)) = aop(6*n); 
thetar(l) = aop(H-6*n); thetar(n 1(1,1)) = aop(7*n); 
for k = 2;(nl-l) 

taut = (2*t(k)-tau)/tau; 
k 

rr(k) = interpl(xt,aop(l:n),taut,'spline'); 
vr(k) = interpl(xt,aop(n+l ;2*n),taut,'spline'); 
gmr(k) = interpl(xt,aop(2*n+l:3*n),taut,'spline'); 
massr(k) = interpl(xt,aop(3*n+l:4*n),taut,'spline’); 
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if massr(k) > massr(k-l) 

massr(k)=massr(k-1); 

end; 

thrustr(k) = interpl(xt,aop(4*n+l:5*n),taut,'spline'); 
if thrustr^) < 0 

thrustr(k) = 0; 

end; 

alphar(k) = interpl(xt,alpha,taut,'splme'); 
thetar(k) = interpl(xt,aop(6*n+l:7*n),taut,'spline’); 

end; 

figure(l) 

erad = x(;,l)'-rr; evel = x(:,2)'-vr;eniass = x(:,4)'-niassr; 

subpIot(2,l,l),plot(t,(x(;,l)'-rr),'b',t,(x(:,2)’-vr),’r’); 

title('Radius Error (Blue) Velocity Error (red)'); 
subplot(2,l,2),plot(t,(x(:,3)’-gmr),'b',t,(x(:,4)'-massr),'r'); 
title('Ganima Error(Blue) Mass Error(Red)'); 

figure(2) 

subplot(2,2,l),plot(t,rr,'r');title('Radius vs Time'); 
subplot(2,2^),plot(t,vr,'r');title('Velocity vs Time'); 
subplot(2^,3),plot(t,gmr,'r*);titleCGamma vs Time’); 
subplot(2,2,4),plot(t,massr,'r');title('Mass vs Time'); 

figure(3) 

subplot(2,l,l)»plot(t5thrustr,'r');tide('Thrust vs Time'); 
subplot(2,l,2),plot(t,alphar,'r');title('Alpha vs Time'); 

band = max(rr)-min(rr); 

rprop = x(:,l); vprop =x(:,2); gmprop = x(:,3);mprop =x(;,4); 
thprop = x(;,5); 

Appendix G - Orbres.m 


global aop costfii 
[costfh,g]=dorbcrit(aop) 

costfii 

r=aop(l:n) 

V = aop(n+l:2*n) 
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gamma = aop(2*n+l:3*n) 
mass = aop(3*n+l;4*n) 

T = aop(4*n+l:5*n) 
alpha = aop(5*n+l:6*n) 
theta = aop(6*n+l:7*n) 
tau = aop(7*n+l) 
y = (tau/2)*(x+l); 

Appendix H - Lobatto.m 

% here are the files for lobatto points. 

% [x w]=lobatto(n) n= number of points. 


function [x,w] = lobatto(n,a,b) 

% [x w] = lobatto(n) or [x w] = lobatto(n,alpha,beta): 

% 

% Computes abscissa and weights for the n-point Gauss-Jacobi-Lobatto 

% quadrature formula using the method of Gene H. Golub, Some modified 

% matrix eigenvalue problems, SIAM Rev. 15 (1973) 318-334. Another early 
% algorithm for this is by David Galant, An implementation of Christoffel's 
% formula in the theory of orthogonal polynomials. Math. Comp. 25 (1971) 

% 111-113. All such Algorithms should be "reviewed", in light of recent 

% improvements in tqr and Cholesky LR algorithms. But, this algorithm 
% "ain't bad". 

% Copyright (c) 23 August 1997 by Bill Gragg. All rights reserved. 

% lobatto calls mxt, mxtj and tqr. 


% begin lobatto 


if nargin < 2 
a = 0; b = 0; 


end 

m = 2'^(a + b + l)*beta(a+l ,b+l); us = a == b; 
n = n - 1; [a b] = mxli(n,a,b); T = mxt(a,b); 

I = eye(n); e = zeros(n, 1); e(n) = 1; 

c = (T + I)\e; c = c(n); d = (T - I)\e; 

d = d(n); e = c-d; c = (c + d)/e; 

d = sqrt(2/e); a(n+l) = c; b(n) = d; 


[xu]=tqr(a,b); u = u’; 


w = m*u.''2; 
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% "Purify" formulas in the ultraspherical case, 
if us 

X = (x - flipud(x))/2; w = (w + flipud(w))/2; 
end 

% end lobatto 

Appendix I - 


function [Dn,x,w]=Diff(n); 

% This function calculates the differentiation matrix Dn that is 
% obtained by differentiating the Legendre Polynomials at the legendre- 
%Lobatto points. It's zero on the main diagonal except at l=k=l, where 
% Dn(l,l)= n(n+l)/4; and at l=k=n; where Dn(n,n)=-n(n+l)/4. 

% n= no of Lobatto points. For the other points 1 (~=)k, we have 
% DnO,k)= Ln(xl)/Ln(xk)*(l/xl-xk). 


[x w]=lobatto(n); 
x=soit(x); 

% initialize Dn 
Dn=zeros(n); 

Dn(l,l)=-(n-l)*n/4; 

Dn(n,n)=n*(n-1 )/4; 

% Calculate the legendre polynomials at xi. 
p=0*eye(n); 

for i=l;n; s=x(i); p(i,l)=l; p(i,2)=s; 

for j=2;n-l; p(ij+l)=((2*j-l)*s*p(ij)-(j-l)*p(ij-l))/j; end; end; 


% Fill out the rest of matrix Dn. 
for l=l:n; for k=l;n; 
ifl-^k, 

Dn(l,k)=p(l,n)/(p(k,n)*(xa)-x(k))); 

end; 

end;end; 

Appendix J - Tfoptm 

% This routine calculates the minimum of the cost function given by 

% the function crit. 

clear; 

tic; 
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global n 

global Dn X w 

global B ve rho aop aref 

n=24; % The number of lobatto points. 

%m0 = 408420; a=2200; 
m0 = 3000;a=500; 

%m0=3000;a=10000; 

isp = 300; gn = 9.81; cd = 2.35; 

rok = 6378+300; rhoO = 1.87*10^-11); rho =1; 

vO = sqrt((3.986005*10^5)/rok)*1000; rO = rok * 1000; 

B = 2 * mO / (r0*rho0*a*cd); 


ve = isp*gn/v0; 


Tguess = 5.0; 

[Dn,x,w]=dififin(n); 

%%% the true solutions, u is the control and z is the state variable. 
y=(x+l)/2; 


aop(l;n)= ones(n,l); 
aop(n+l:2*n) = ones(n,l); 
aop(2*n+l:3*n) = l.l*ones(n,l); 
aop(3*n+l:4*n) = ones(n,l); 
aop(4*n+l:5*n) = Tguess*ones(n,l); 
aop(5*n+l:6*n) = 0*ones(n,l); 
aop(6*n+l:7*n) = zeros(n,l); 
aop(7*n+l) =112; 
taup=aop(7*n+l); 
options( 13 )=5 *n+7; 
options(14) = 400*(7*n+l); 
options(4)= 10^(-4);%5 
options(3)= 10^(-5);%5 
options(2)= 10^(-5);%5 
options(l)=l; 

%options(18)=l; 
options(16) = 10'^(-6);%7 
vub=[ ]; 
vlb=[ ]; 
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aop=constr(tfcrit',aop,options,vlb,vub); 

% a is for the state and b is for the control. 

rp=aop(l:n); 

vp = aop(n+l:2*n); 

gammap = aop(2*n+l;3*n); 

massp = aop(3*n+l:4*n); 

Tp = aop(4*n+l;5*n); 
alphap = aop(5*n+l;6*n); 
thetap = aop(6*n+l:7*n); 
tau = aop(7*n+l); 

toe; 

tic; 

orbres; 

orbeonv; 

Appendix K -Tfcritm 

function [costfii,g]=tfcrit(aop); 

%% This function calculates the cost function that is to be minimized 
%% and the state constraints. 

global nx w Dn; 
global B ve rho aref; 

% global costfii; 

% Set up the cost function. 
fori=l:n 
fii(i)=aop(4*n+i); 
end; 

costfii= l/(2)*sum(w.*fa'); 
tau = aop(7*n+l); 


% Set up the state constraints. 
fori=l;n 

rt = aop(i)*(6378+300); 
ifrt>=6728 

rref = 6778; h=58.2; 
rhoref = 2.62e-12; 
elseif rt >= 6628 

rref= 6678; h =50.3; 
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rhoref= 1.87e-ll; 
else rref = 6528; h= 37.5; 
rhoref=2.41e-10; 

end; 

rhon = rhoref*exp(-l*(rt-iTeQ/h); 
rho = rhon/1.87e-ll; 

G(i) = l/(aop(i)^2); 

D(i) = rho*(aop(n+i)^2); 

Art = (aop(4*n+i)*cos(aop(5*n+i))-D(i))/(aop(3*n+i)*B); 
Ant = aop(4*n+i)*sin(aop(5*n+i))/(aop(3*n+i)*B); 

V(i) = ve; 


% Theta 

g(i+4*n) = (2/tau)*sum(Dn(i,:).*aop(6*n+l;7*n))-(aop(n+i)/aop(i)*cos(aop(2*n+i))); 
% Mass 

g(i+3 *n)= (2/tau)* sum(Dn(i,;). *aop(3 *n+1:4 *n))+(aop(4 *n+i)/(ve*B)); 

% Gamma 

g(i+2*n)= (2/tau)*sum(Dn(i,;).*aop(2*n+l :3*n)); 

g(i+2*n^ g(i+2*n)-(((aop(n+i)^2)/aop(i))-G(i))*(cos(aop(2*n+i))/aop(i+ii)); 
g(i+2*n^ g(i+2*n)- Ant/aop(n+i); 

% Radius 

g(i)= (2/tau)*sum(Dn(i,:).*aop(l :n))-aop(n+i)*sin(aop(2*n+i)); 

% Velocity 

g(i+n)= (2/tau)*sum(Dn(i,:).*aop(n+l ;2*n))+(G(i)*sin(aop(2*n+i)))-Art; 


g(i+5*n+7) = - aop(4*n+i); 
g(8*n+7+i) = aop(4*n+i)-5.0; 

% g(i+9*n+7) = l-aop(i); 
%g(i+7*n+7) = aop(i)-l; 
g(i+6*n+7) == aop(5*n+i)-pi; 
g(i+7*n+7) = -pi-aop(5*n-i-i); 
end; 

g(5*n+l) = aop(l)-aop(n); 
g(5*n+2) = aop(2*n+l)-aop(3*n); 
g(5*n+3) = aop(6*n+l); 
g(5*n+4) = aop(n+l)-aop(2*n); 
g(5*n+5) = aop(3*n+l)-l; 
g(5*n+6) = aop(l)-l; 
g(5*n+7) = aop(7*n+l)-l 12; 

g=g'; 
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Appendix L - Notoptm 


% This routine calculates the minimum of the cost function given by 

% the function crit. 

clear; 

tic; 

global n 

global Dn X w 

global B ve rho aop aref 

n=24; % The number of lobatto points. 

%m0 = 408420; a=2200; 
mO = 3000; a =500; 

%m0=3000;a=10000; 

isp = 300; gn = 9.81; cd = 2.35; 

rok = 6378+300; rhoO = 1.87*10^-11); rho =1; 

vO = sqrt((3.986005*10^5)/rok)*1000; rO = rok * 1000; 

B = 2 * mO / (r0*rho0*a*cd); 


ve = isp*gn/v0; 


Tguess = 25.0; 
[Dn,x,w]=diflfin(n); 


aop(l:n)= ones(n,l); 
aop(n+l;2*n) = ones(n,l); 
aop(2*n+l;3*n) = l.l*ones(n,l); 
aop(3*n+l:4*n) = ones(n,l); 
aop(4*n+l:5*n) = Tguess*ones(n,l); 
aop(5*n+l;6*n) = 0*ones(n,l); 


taup=112.6; 
options(l 3)=4*n+5; 
options(14) = 400*(7*n+l); 
options(4)=l 0^-6);%5 
options(3^10^(-8);%5 
options(2^10^(-8);%5 
options(l)=l; 
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%options(18)=l; 
options(16) = 10^(-7);%7 
vub=[]; 
vlb=[ ]; 

aop=constr('notcrit',aop,options,vlb,vub); 

% a is for the state and b is for the control. 

rp=aop(l;n); 

vp = aop(n+l;2*n); 

gammap = aop(2*n+l:3*n); 

massp = aop(3*n+l:4*n); 

Tp = aop(4*n+l:5*n); 
alphap = aop(5*n+l;6*n); 

tau= 112.6; 

toe; 

tic; 

notres; 

notconv; 

toe; 

Appendix M - Notcritm 

function [costfii,g]=notcrit(aop); 

%% This function calculates the cost function that is to be minimized 
%% and the state constraints. 

global nx wDn; 
global B ve rho aref; 

% global costfii; 

% Set up the cost function, 
for i=l:n 

fii(i)=aop(4*n+i); 

end; 

costfii= l/(2)*sum(w.*fii'); 
tau = 112.6; 


% Set up the state constraints. 
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fori=l:n 

rho = density(aop(i)); 


G(i) = l/(aop(i)^2); 

D(i) = rho*(aop(n+i)^2); 

Ast = (aop(4*n+i)*cos(aop(5*n+i))-D(i))/(aop(3*n+i)*B); 
Ant = aop(4*n+i)*sin(aop(5*n+i))/(aop(3*n+i)*B); 

VCi) = ve; 


% Theta 
% Mass 

g(i+3 *n)= (2/tau)*sum(E)n(i,;). *aop(3*n+l ;4*n))+(aop(4*n+i)/(ve*B)); 

% Gamma 

g(i+2*n)= (2/tau)*sum(Dn(i,:). *aop(2*n+l :3 *n)); 
g(i+2*n)=g(i+2*n)-(((aop(n+i)^2)/aop(i))-G(i))*(cos(aop(2*n+i))/aop(i+n)); 
g(i+2*n)= g(i+2*n)- Ant/aop(n-t-i); 

% Radius 

g(i)= (2/tau)*sum(Dn(i,;).*aop(l;n))-aop(n+i)*sm(aop(2*n+i)); 

% Velocity 

g(i+n)= (2/tau)*sum(Dn(i,:). *aop(n+l :2*n))+(G(i)*sin(aop(2*n+i)))-Ast; 


g(i+4*n+5) = - aop(4*n-H); 
g(5*n+5+i) = aop(4*n+i)-5.0; 
g(i+6*n+5) = aop(5*n+i)-pi; 
g(i+7*n+5) = -pi-aop(5*n+i); 
end; 

g(4*n+l) = aop(l)-aop(n); 
g(4*n+2) = aop(2*n+l)-aop(3*n); 
g(4*n+3) = aop(n+l)-aop(2*n); 
g(4*n+4) = aop(3*n+l)-l; 
g(4*n+5) = aop(l)-l; 

g=g'; 

Appendix N - Noaoptm 

% This routine calculates the minimum of the cost function given by 

% the function crit. 

clear; 

tic; 
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global n rhom vm 
global Dn X w 

global B ve rho aop aref band 
global tg ag taup 

n=24; % The number of lobatto points, 
band = 0.01; 

rhom = density(l+band/2); 
vm= sqrt(l/(l+band/2)); 

m0 = 408420; a=2200; 
mO = 3000; a =500; 

%m0=3000;a=10000; 

isp = 300; gn = 9.8067; cd = 2.35; 

rok = 6378+300; rhoO = 1.87*10^(-11); rho =1; 

vO = sqrt((3.986005*10^5)/rok)*1000; rO = rok * 1000; 

B = 2 * mO / (rO*rhoO*a*cd); 


ve = isp*gn/v0; 


Tguess = 20; 

[Dn,x,w]=difi&n(n); 

aop(l:n)= ones(n,l); 

aop(n+l:2*n) = ones(n,l); 

aop(2*n+l:3*n) = .5*ones(n,l); 

aop(3*n+l:4*n) = ones(n,l); 

aop(4*n+l:5*n) = Tguess*ones(n,l), 

%aop(4*n+8:5*n) = ones(length(aop(4*n+8;5*n)),l); 

taup=112.6; 

options(13)=4*n+5; 

options(14) = 400*(7*n+l); 

options(4)=l(r(-7);%5 

options(3)=l 0''(-8);%5 

options(2)=10^-9);%5 

options(l)=l; 

%options(18)=l; 
options(16) = 10^-8);%7 
vub=[ ]; 
vlb=[]; 

aop=constr('noacrit',aop,options,vlb,vub); 
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% a is for the state and b is for the control. 

rp=aop(l:n); 

vp = aop(n+l;2*n); 

ganimap = aop(2*n+l;3*n); 

massp = aop(3*n+l:4*n); 

Tp = aop(4*n+l:5*n); 


tau = taup; 

toe; 

tic; 

noares; 

noaconv; 

toe; 

Appendix O - Noacritm 

function [costfii,g]=noacrit(aop); 

%% This function calculates the cost function that is to be minimized 
%% and the state constraints. 

global nx wDnband; 
global B ve rho aref; 
global rhom vm taup; 

% global costfii; 

% Set up the cost function. 
i=l:n; 

% fii=aop(4*n+i)/(rhom*vm^2); 

fii=aop(4*n+i); 
costfo= l/(2)*sum(w.*fii'); 
tau= 112.6; 


% Set up the state constraints, 
for i=l:n 


rho = density(aop(i)); 
G(i) = l/(aop(ir2); 
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D(i) = rho*(aop(n+i)^2); 

Ast = (aop(4*n+i)*l-D(i))/(aop(3*n+i)*B); 
Ant = aop(4*n+i)*0/(aop(3*n+i)*B); 

V(i) = ve; 


% Theta 
% Mass 

g(i+3 *n)= (2/tau)*suni(Dn(i,:). *aop(3 *n+l :4*n))+(aop(4*n+i)/(ve*B)); 

% Gamma 

g(i+2*n)= (2/tau)*sum(Dn(i,:).*aop(2*n+l;3*n)); 
g(i+2*n)=g(i+2*n)-(((aop(n+i)^2)/aop(i))-G(i))*(cos(aop(2*n+i))/aop(i+n)); 
g(i+2*n)= g(i+2*n)- Ant/aop(n+i); 

% Radius 

g(i)= (2/tau)*sum(Dn(i,:).*aop(l :n))-aop(n+i)*sin(aop(2*n+i)); 

% Velocity 

g(i+n)= (2/tau)*sum(Dn(i,:).*aop(n+l:2*n))+(G(i)*sin(aop(2*n+i)))-Ast; 


g(i+4*n+5) = - aop(4*n+i); 
g(5*n+5+i) = aop(4*n+i)-20.0; 
g(i+6*n+5) = aop(i)-(l+band); 

end; 

g(4*n+l) = aop(l)-aop(n); 
g(4*n+2) = aop(2*n+l)-aop(3*n); 
g(4*n+3) = aop(n+l)-aop(2*n); 
g(4*n+4) = aop(3*n+l)-l; 
g(4*n+5) = aop(l)-l; 


r=g'; 


Appendix P - Tqr.m 

% Total flops (scalar case, see csgn); TBC 


% Problem. 
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% 1. Compare this function experimentally with csgn. Compare with regard 

% to both execution time and numerical stability. Is matlab computing 
% sign correctly? 


function [lam,U] = tqr(a,b,U) 

% [lam u] = tqr(a,b) or [lam U] = tqr(a,b,U); 

% 

% [lam u] = tqr(a,b): 

% 

% The column lam contains the eigenvalues of the Hermitian tridiagonal 
% matrix T = mxt(a,b) computed by one version of the (real symmetric) tqr 
% algorithm with Wilkinson's shift. The column u contains the first 
% elements of the eigenvectors of T normalized to be nonnegative and such 

% that the eigenvectors are unit vectors. In practice this is an 0(n^2) 

% process. Ifu is omitted only the eigenvalues are computed. The 
% computed eigenvalues are real and are sorted to be nonincreasing. 

% 

% [lam U] = tqr(a,b,U); 

% 

% This replaces the input U by UV with V a matrix of orthonormal eigen- 
% vectors of T. Ifthe input U is I the output U is V. Ifthe input U is 

% unitary with AU = UT then the output U is unitary with AU = UD and D = 

% diag(lam). 

% 

% If the input U is e(l)' the output U is u'. If the input U is 

% [e(l)'; e(n)'] the output U is [u'; v'] with v the column of last 

% elements of the normalized eigenvectors. Ifthe subdiagonal elements of 
% T are all nonzero then the elements of v alternate in sign, at least 
% mathematically. 

% Copyright (c) 2 February 1991 by Bill Gragg. All rights reserved. 

% Revised 15 July 1994. 

% tqr calls sgn. 


% begin tqr 

% Ensure that T is Hermitian and shift b down one unit. 
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a = rea](a); n = length(a); b = [0;b(;)]; b = b(l;n); 


% Initialize U if required and execute a diagonal unitary similarity 

% transformation to make T have nonnegative subdiagonal elements. 

if nargout > 1 

if nargin < 3 

U = zeros(l,n); U(l) = 1; 
end 

u = sgn(b); u = cumprod(u); U = U*diag(u); 
end 

b = abs(b); 

% Scale the matrix up by a power of two to give nearly the widest 
% possible exponent range. 

scale = norm([a; b*sqrt(2)]); scale = 2'^(1024 - ceil(log2(scale))); 
a = a*scale; b = b*scale; 


format compact % Temporary statements 

ibaxscale = max(abs([a; b])); % for display. 

% "Dotqr". 


form = n:-l;l 


% disp(m) % Temporary display statement. 

% Compute the mth eigenvalue. 

for its = 0; 10*n % its is the iteration index. 

% Split the matrix if possible. This is also the termination 

% test. 


fork = m:-l:l 


ifk> 1 

tol = abs(a(k-l)) + abs(a(k)); 
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if tol + b(k) = tol 
b(k) = 0; break 
end 


end 

end 


ifk = m 


break %b(m) = 0. a(m) is an eigenvalue, 
else 


if its = 10*n 

error(tqr iteration did not taminate in lOn steps!') 
end 

% Compute Wilkinson's shift w as a perturbation of the 

% Rayleigh shift r = a(m). As the algorithm converges 

% c = b(m) —> 0. 

r = a(m); c = b(m); d = (r - a(m-l))/2; s = abs(d); 

ifc<s 

s = c/s; t = 1 + sqrt(l + s*s); t = c*s/t; % t < c; 
else 

s = s/c; t = s + sqrt(l + s*s); t = c/t; % t < c; 
end 

ifd>0 
w = r +1; 
else 

w = r -1; 
end 

% Take a step of the tqr algorithm. There are many ways to 

% implement the inner loop. We recently found the fastest 

% known stable form in terms of flops. The form given here 

% is elegant. 


c=l; s = 0; p = w-a(k); t = p; 
forj = k:m-l 
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% Compute the two by two reflector stably and update b(j). 

oldc = c; oldt = t; q = b(j+l); u = abs(p); 
ifq<u 

V = q/u; r = sqrt(l + v*v); b(j) = u*r*s; 

u = sgn(p); c = u/r; s = v/r; 

else 

V = p/q; r = sqrt(l + v*v); b(j) = q*r*s; 

u=l; c = v/r; s = u/r; 

end 

% Update p, t, a(j), and U(;,j:j+1) if required. 

p = c*(w - aO’+l)) - s*q*oldc; t = c*p; 
a(j) = a(j+l) +1 - oldt; 

if nargout > 1 

i=j:j+l; U(:,i) = U(;,ir[-C s; s c]; 
end 

end 

% Update b(m), a(m), and U(:,m) if required. 

b(m) = abs(p)*s; a(m) = w -1; c = sgn(p); 

% disp(b(m)/maxscale) % Temporary display statement. 

if nargout > 1 
U(:,m) = -U(:,m)*c; 
end 

end 

end 

% pause(3) % Temporary pause statement, 
end 

% Sort and prepare the output. 

[a p] = sort(-a); lam = - a/scale; 
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if nargout > 1 

U = U(:,p); u = U(l,:)’; 

if nargin < 3 
u = abs(u); U = u'; 
else 

u = sgn(u); U = U*diag(u'); 
end 

end 

% end tqr 

Appendix Q - Mxtm 

% Problems. 

% 1. Relate T = nixt(a,b), with [ab] = inx^(n,l/2X with the negative second 
% dififCTence matrix S = mxt(c,d), with [c d] = mxs(n). 

function T = mxt(a,b,c) 

% T = mxt(a,b,c) or T = mxt(a,b): 

% 

% T = mxt(a,b,c) is the TRIDIAGONAL MATRIX with diagonal elements a(l :n), 
% subdiagonal elements b(l :n-l) and superdiagonal elements c(l :n-l), 

% 

% T = mxt(a,b) is the HERMITIAN tridiagonal matrix with diagonal elements 
% real(a(l :n)) and subdiagonal elements b(l :n-l). 

% Copyright (c) 1 December 1990 by Bill Gragg. All rights reserved. 

% Revised 21 November 1992. 

% mxt calls no extrinsic functions. 

% begin mxt 

if nar^ < 3 
a = real(a); c = b'; 
end 

n = length(a); b=b(l:n-l); c = c(l:n-l); z = zeros(n-l,l); 
ifn<500 


B = diag(b); B = [z'0;Bz]; C = diag(c); C = [zC;0z']; 
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T = diag(a); T = T + B + C; 
else 

T = zeros(n); 
fork= l:n-l 

T(k,k) = a(k); T(k+l,k) = b(k); T(k,k+1) = c(k); 
end 

T(n,n) = a(n); 
end 

% end mxt 

Appendix R - Mxtj.m 


function [a,b] = nixtj(n,alpha,beta) 

% [a b] = nixtj(n,alpha,beta), [a b] = nixQ'(n,alpha), [a b] = nixq'(n), 

% T = nix^(n,alpha,beta), T = nixtj(n,alpha) or T = nixtj(n): 

% 

% nix^(n,alpha,beta); T = nixt(a,b) is the Jacobi matrix whose characteristic 
% polynomial p is (a nonzero scalar multiple of) the nth JACOBI polynomial. 

% The eigenvalues of T are the abscissas of the nth order Gauss-Christoffel 
% quadrature formula for the weight function ((1 - t)^alpha)((l + t)^beta) on 
% the interval - 1 < t < 1. The Gauss-ChristofFel weights are m(0) times the 
% squares of the first elements of the normalized eigenvectors of T, where 
% m(0) = b(0)^2 = B(alpha + l,beta + l)2^alpha + beta -1) is the total mass. 

% B is the beta function. The weight function is positive and integrable if 
% alpha + 1 > 0 and beta + 1 > 0. 

% 

% mx^(n,alpha) takes beta = alpha, p is the nth ULTRASPHERICAL polynomial, 
% with weight function (1 - t‘^2)^alpha on the interval -1 < t < 1. Special 
% cases are the CHEBYSHEV polynomial of the FIRST KIND, with alpha = -1/2, 
% and of the SECOND KIND, with alpha = 1/2. 

% 

% mx^(n) takes alpha = beta = 0. p is the nth LEGENDRE polynomial, with 
% weight function w(t) = 1 on the interval - 1 < t < 1. The quadrature 
% formula here is originally due to Gauss. Christoflfel generalized Gauss' 

% formula to a wide class of weight functions. Because of this the Gauss- 
% Christoffel weights are usually called Christoffel numbers. 
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% Copyright (c) 2 Febmaiy 1991 by Bill Gragg. All rights reserved. 

% mxtj calls mxt. 

% begin mx^ 

ifnargin<2 alpha = 0; end; ifnargin<3 beta = alpha; end 
a = alpha; b = beta; c = a + b; d = b-a; 
s(l) = d/(c + 2); t(l) = (a + l)*(b + l)/(c + 2)^2/(c + 3); 
ifn>2 
d = c*d; 

n = (2:n)'; m = 2*n; mm = ni-l; mp = m+l; 
s(n) = d./(c + m)./(c + (m - 2)); 

t(n) = n.*(a + n).*(b + n).*(c + n)./(c + mm)./((c + m).'^2)./(c + mp); 
end 

a = s(:); b = 2*sqrt(t(:)); 
ifnargout<2 a = inxt(a,b); end 
% endmx^ 
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